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MATHEMATICAL AND COMPUTATIONAL STUDIES OF THE STABILITY 
OF AXISYMMETRIC ANNULAR CAPILLARY FREE SURFACES 

by Norman Albright, Paul Concus, and Ilkka Karasalo 


We present here the results of our mathematical and computational 
studies of the stability of a liquid in a rotationally symmetric 
container subject to gravitational and surface forces. Of specific 
interest is the case for which the' contact angle is zero, or nearly 
zero. The application of primary concern is that of stability in a 
vertical right circiilar cyclindrical container with a concave 
spheroidal bottom, for the case in which the volimie of liquid is 
sufficiently small so that liquid lies only in an annular region of 
the container. Numerical computations are presented for a container 
used for the storage of liquid fuels in National Aeronautics and Space 
Administration Centaur space vehicles, for which the axial ratio of 
the container bottom is 0,724. 

Our studies consist of several self-contained parts, which are 
discussed independently in the appendices that follow. • In Appendices I 
and II are derived the mathematical results on which our static-analysis 
computations are based. These results are concerned with the conditions 
for the contained liquid to be in stable equilibriinn. Of particular 
interest is the case of zero contact angle, which has pre'viously not 
received adequate mathematical treatment in the literature. 

A" configuration is in stable equilibrium if and only if it 
strictly minimizes the sum of the surface and gravitational potential 
energies, among all nearby configurations ■with the same liquid volume. 
With a suitable choice of variables, the problem of • stability may be 
approached as a variable-endpoint problem in the calculus of variations. 
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Using this approach conditions are obtained that distinguish 
between stable and unstable cases in. a satisfactory way-j 
if. the contact angle is not zero. If the contact angle is zero, the 
formal limits of the endpoint conditions depend crucially on whether 
or not the curvatures of the equilibrium interface and the container 
wall coincide at the three phase contact lines. If they do not, the 
limiting endpoint conditions will be of the fixed type. 

It is shown that for zero contact angle the stability 
criteria based on the fixed end-point conditions apply if only the 
analytic continuation of the equilibrium liquid-vapor interface does 
not penetrate the container walls at the three-phase contact lines. 

If it does penetrate, then the configuration will be unstable regardless 
of the conditions on the second variation of the total potential 
energy. 

A computational study based on the mathematical results in Appendices 
I and II is carried out in Appendix III. The critical Bond number for 
stability is calculated for the Centaur space vehicle tank as a 
function of the liquid volume for contact angles of 0®, 1°, 2°, and 4®. 
For the zero-contact-angle case, the specially derived end-point condi- 
tions are used, while for the latter three cases the usual variable 
end-point conditions are used . The zero-contact angle results are 
found to be consistent with the limiting ones for the nonzero angles. 
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Of particular interest for the case of zero contact angle is the 
appearance, for this problem, of a critical liquid volume marking a 
qualitative change in the nature of the solution. For liqiiid voliimes 
less than the critical one the stability limit is determined by the 
fixed end-point conditions, whereas for liquid volumes larger than the 
critical one the stability limit is determined by the non-penetration 
criteria discussed above. This feature corresponds to the existence 
of (unstable) equilibrium coiiflgurations for Bond nunibers larger in 
magnitude than the critical one only for the small volume case. The 
critical liqttid volume for the' Centaur tank corresponds to a fill 
hei^t of 0.5031 a, where a is the radius of the tank. 

The numerical results presented in Appendix III have been found 
to be consistent with preliminary e3q>erimental results obtained by 
E. F. Symons at the NASA Lewis Research Center. 

In Appendix IV small-amplitude^periodic sloshing modes are 
calculated for the container configuration studied in Appendix III, 
for zero contact angle. As must be the case, for a conservative 
mechanical system, the critical Bond numbers for stability obtained 
from the dynamic analysis are found to agree with those calculated 
from the static analysis. Agreement is also found on the value of 
the critical liquid volume for which equilibrium configurations 
can exist nearby the critically stable one. 

Oscillation frequencies or growth rates are calculated for 
several Bond numbers and liquid volumes, for normal modes having up 
to six angular nodes. The computations indicate that for liquid 
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volumes smaller than the critical one, the sloshing modes without 
radial nodes generally become unstable before those with one or more 
radial nodes, as the magnitude of the Bond number is increased. The 
calculations for larger liquid volumes indicate that all modes become 
unstable together, in agreement with the cessation of existence of 
a nearby equilibrium for this case. 

The above results are depicted graphically, and growth rates for 
the dominant unstable mode for a liquid volume corresponding to a 
dimensionless mean fill height of 0,30 are given for three different 
liquids. 
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ABSTRACT 

The stability, in terms of minimum static energy, of axisymmetric annular 
liquid interfaces in axisymmetric containers is studied mathematically. The 
sufficiency of stability conditions based on the Jacobi accessory minimization 
problem with respect to volume conserving, "weak" perturbations is proved in 
two cases: Firstly, sufficient conditions' for stability in the general case 
of nonzero liquid-wall contact angle are considered. Secondly, the formal 
limit of these conditions, as the contact angle tends to zero is proved to be 
correct, provided that the curvatures of the unperturbed liquid-vapor inter- 
face and the container wall are not equal at their contact lines. 

^^^^(P*'6ceding page blank 
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1. FORMULATION OF THE PROBLEM 

In this paper we study mathematically the stability of certain configur- 
ations of liquid contained in an axially S3nnmetric tank in a .gravitational 
field directed along the axis of symmetry. We require that the tank shape 
and the liquid volume are such that the liquid- vapor interface is annular, 
i.e. it does not intersect the axis of S3nmnetry, cf., fig. 1: 


Container 

wall 


XBL 766-3054 

Figure 1: Example of permissible liquid-tank configuration 

and the associated cylindrical coordinate system. 

A configuration is one of stable equilibrium if and only if it 




9 


strictly minimizes the total static potential energy of the system, given 
by 


E = a(A. - cos Y A ) + E , 
f w' g ’ 


( 1 . 1 ) 


compared' to all nearby configurations with the same liquid volume V. Here 

a ^ 0 and 0 < Y < IT are constants determined by the physical properties of the 

liquid and the wall, and are the areas of the liquid- vapor and liquid- 

container interfaces, respectively, and E is the gravitational potential 

6 

energy of the liquid. 

We shall study below the behavior of the energy (1.1) when the liquid is 
perturbed slightly (but not necessarily rotational symmetrically) from a 
configuration of rotational symmetry. For our purposes, a parametric arc- 
length, normal displacement representation of the surfaces (cf . , e.g. Reynolds, 
Saad, Satterlee [ 2 ]) offers some advantage. Hence we let the unperturbed 
surface be described by the equations 


r = R(s) ; 0 < 4) < 2ir- 
z = Z (s) ; < s < s^ 


( 1 . 2 ) 


in the polar co-ordinate system of figure 1 where the' parameter s is' the* arc- 
length along the curve of intersection between the liquid-vapor surface and any 
plane (j> = constant. Then the equations 


r = R(s) - n(s,<J>)Z’ (s) ; 0 < (f) < 2 tt 
z = Z (s) + q(s,4))R’ (s) ; s^((j)) < s < s^(^) 


(1.3) 


describe a surface obtained by moving each point of the surface (1.2) the 
distance r|(s,({)) in the direction of the principal normal. In general, since 
we want the perturbed surface (1.3) to intersect the container walls, the 
functions R and Z of (1.2) must be continued in the s-direct±on to some open 
interval containing [ s^ , s^ 3 , cf. figure 2. A convenient way of doing this, 
which we will use in the sequel, is provided by the differential equations 
(1.11) below. 
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Figure 2: The arc— length, normal displacement coordinates; 

Intersection between the surfaces and some plane 
(p - constant. 


Similarly, we let the container wall in some neighborhoods of the unper- 
turbed contact lines (s = s and s = s. in (1.2)) be given by 

o 1 


r = E,(s) - w(s)Z’(s) 
z = Z(s) + w(s)R’ (s) 


0 < <|) < 2tt . 


(1.4) 


((1.4) is not adequate if the contact angle, y> between the unperturbed surface 
and the wall is tt/ 2. To avoid unnecessary detail we therefore assume y ^ ir/2 
in the following). 

For any equilibrium configurations , the functions R and Z will be suffi- 
ciently smooth for the perturbed surface (1.3) not to intersect itself when 
|ri(s,(j))|, |ti^(s, 4))] and |r|^(s,«j>)l are small. (For this, it is sufficient that 
R"(s) and Z"(s) are bounded in some open interval containing [s^,s^], cf. 
Concus, Crane, Satterlee [5], p. 4-6, while in fact it follows from the 
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differential equations (1.11) below that R and Z may be, continued analytically 
to some open region containing in the complex s-plane). Then the 

increments of the quantities appearing in (1.1) and of the volume may be 
calculated in a straightforward way. One obtains 


<SE(q) = a • (5A_ - cos Y * 6A ) + 6E 
f ' w g 

2Tr s- (^) 


.... 

= I I |f^(Tl,s) - f^(£,s) + fg(jl,s)| dsd<{) 

J 0 -^s^(^) 


-I I |cos Y ■ f^(w»s) - f^(£,s) + fg(w^s)|' dsd(J) 

•' 0 *' A(() 


(1.5) 


■f.r 


6v(n)= 


0 S^(<f>) 


f.^(q,s) dsd<{) 


27T 

-n 


f ^ (w, s ) dsd<j) 


0 


where we have put 


( 1 . 6 ) 


T 

]i=ii(s,(j>) = (Ti(s,<|)), Tig(s,(i)), n^(s,4»)) 
w = w(s) = (w(s) ,w* (s) ,0)^ 

A(j) = the interval (s^(4>) ,s^)u(s^,Sj^((f))) 
and, denoting R = R(s), Z = Z(s), 

f^(n,s) = a ♦ {(R - riz')^(TiJ + (1 + n(R"z’ - z"r’))^) 
+ n?(i + Ti(R"z’ - z"r’))^P^^ 
f^(n,s) = nji + f (R"z* - z"R')|(R - §z') 
fg(H>s) = pg • (z + Y rO • fy(n»s) 


(1.7) 


( 1 . 8 ) 
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where p, the liquid density, is assumed to be constant and g is the gravi- 
tation constant with g > 0 if the gravitation force is directed towards the 
negative z -axis . We note that when R > 0 all of the f imctions f^, f ^ and 

f are smooth (in fact analytic) in the arguments n, ri and rt, in a neigh- 
g - _ s (p- - 

borhood of Ti = ^. 

A first necessary condition for T)' = 0 to be a local minimizer of 6E (We 

shall need to make the meaning of "local" more precise below. For the moment 
we may consider all ri which are continuous and smooth in some open region 
containing s^((|)) < s < s^(<j)), 0 < 4* < 2 tt.) when 6V = 0 is that, for 
some constant X all first-order ri-terms in 6E-X6V vanish. Using the notation 
in (1.5) - (1.8) this necessary condition is expressed by the Euler-Lagrange 
equation 



(^,s) 


3f 

+ 3 ^( 0 ,=) 





= 0 




in s < s 
o - 


< V 


with the additional condition #• 


(1.9) 


cos Y ^ 


i=0,l. J 


Putting B = pg/d, H = X/o, using (1.8) and the identity 
o o 

R'(s)^ + Z’(s)^ = 1 


( 1 . 10 ) 


(1.9) becomes: 


Z" = R’(B Z - H - ZVR) 
o o 

R" = -Z’(B Z - H - Z'/R) 
o o 

with the boundary condition, that the curve 

(R(s),Z(s)) intersects the container walls 

at s = s and s = s^ , under the contact 
o 1’ 

angle y. 


( 1 . 11 ) 
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(without loss of generality we have excluded the contact, angles 7 T - 7 also 
allowed by (1.9)). 

Assuming now that R and Z are such that (1.11) and, consequently, (1.9) 
are satisfied, a second necessary condition will -be that the terms of .second 
order in Ti, and, ri^ will give a non-negative contribution to 6E for all 
n satisfying the volxme constraint. When y > 0, 'we may formulate this condition 
as follows: 

2TT s^ 

•^0 s 

o 


2-jr 

+ I |a^n(s^,(}))^ + a^Ti(s^,(l))^ d(|)|> 0 

n 


0 


for all n such that 


Zir s., 

// 

*^0 •'s 


3f 


V 


3n 


(£,s) * Tl(s,(|)) dsd(|) = 0 


(1.1-2) 


Here we have put 


_ (-1)^ d 


2 

tan Y 


|^|cos Y f^(w(s),s) - f^(j0,s) + fg(w(s),s) - Xf^(w(s),s)| 


s=s . 
1 


i=0,l. (1.13) 


Using (1.7) - (1.8), the second derivatives appearing in the first integral 
in (1.12) are seen to be diagonal 3X3 matrices, and (1.12) reduces to 
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for all ri such that 



A(s)rt(s,<{») dsdcj) == 0 ' 


J 


(1.14) 


where A(s) = ctR(s),- B(s) = a/R(s) and G(s) is smooth (analytic) in [s^,s^], 
see (5.1) below. 

By standard results for symmetric, semibounded quadratic forms in .Hilbert 
space (see e.g. Kato [l] p. 322 and pp. 352-353), (1.14) may be analyzed in 
terms of the eigenvalues and eigenfunctions of an associated self adjoint 
differential generator: 


Ttl - - |j(A(s)H^) - + C(s)Tl 




ins <s<s_,0<d)< 2ir, with the 
boundary conditions, that H should be 
periodic in (f) with period 2ir and 


A(s )n (s ,<|)) = (-l)^a n(s ,(})) , 0 < ^ < 27T, 

J- o i. XI “ ~ 


(1.15) 


i=0,l 


T has a complete, orthogonal system of eigenfunctions of the form 

cos sin 

with associated eigenvalues (in increasing order in the index i) 


l^ifc}i==l,lc=0 


(1.16) 
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determined from the boundary-value problems 


- + C(s)Sy,j^(s) - 

(k=0,l,2.. . 

A(s.)u! (s.) = (-l)^a.y., (s ) ; j=0,l U=l,2,'3... 

j J j‘ J ‘ *17) 

We notice that all eigenfunctions (1.16) but those with k = 0 satisfy 

the side condition in (1.14). Furthermore, since B(s) > ,0 in s < s < s •, 

o — - i 

we see from (1.17) that the eigenvalues are increasing functions of k. 

It then follows that (1.14) holds for all rt in the class of continuous 
functions in < s < s^, 0 < tfi < 2Tr which are periodic in (j) with period 
2 it and have square integrable first derivatives (see e.g. Kato [1], p. 323, 
Cor. 2.3) if and only if 



where, denoting (f 


Xq “ min|g^X^P + 32^20’ ^11 } - 
Zn s. 

2 2 

f (s,(j>)g(s,^)dsd(j), and $2 sre the solutions 

: 


3^(Uio>^)o ^2^^20’^^o “ ® 


%^^10’*^10^o ^2^’^20’’^20^o ^ 


(1.19) 


(with 32 ~ 0 if solutions are non-unique) . 

The purpose of the present paper is to provide a theoretical complement 
to a computational study of this kind (Concus, Karasalo [7] ) . There the 
necessary condition (1.18) is used to distinguish between stable and' unstable 
equilibrium surfaces in the limiting case Y = 0. We will therefore here 
first show, that when Y ^ 0? (1.18) with strict inequality is. a sufficient 
condition for our constrained minimization problem, when T] is allowed to 
vary in a "weak” type of neighborhood (cf. e.g. Bolza [3], p. 68-70) of 
ri = 0. Then, as the main result (which may be new cf. Gillette [4], p. 23), 
we will show that (1.18) with strict inequality and the eigenvalues computed 
from (1.17) but the boundary condition 
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- 0 ; j - 0,1 : u.ao) 

is a sufficient condition for the constrained minimization problem (again, 

■in a "weak" sense [3] ) when y - 0 provided that the curvatures of the unper- 
turbed surface and the container wall are not equal at the unperturbed contact 
lines . (This additional condition on the unperturbed surface in fact turns 
out to be necessary (see (5.3) below) in order for^ the boundary conditions 
of (1.17) to converge formally to (1,20) as y 0*) 
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2. SOME NOTATION AND LEMMAS 

We denote by S the closed domain s^(4>) < s < s^((j)), 0 < (J) < 2 it in the 
(s,<f))-plane. We require that a permissible perturbing function 
q= ri(s,(})) in (1.3) should have the following regularity properties: 

a) ri is continuous in S and periodic in (j) with period 277. 

b) T]^ and are continuous in S except possibly at finitely many points 

or finitely many piecewise smooth curves with finite length. In 

particular, ri is piecewise continuous as function of s for all (}>. 
s 

c) d(ri) “ supin(s,(j))| + sup]rig(s,0)l + sup (ri^(s,(j)) | < “ (2.1) 

where the supremum is taken over all points of S where and are 
continuous . 

d) The functions s^ ((()); i=0,l are continuous, and such that s^((J)) - s^, i=0,l, 

change sign at most finitely many times in 0 <‘ <|) < 277. s^'.((j)); i=0,l, are 

continuous except possibly at finitely many points in 0 < (j) < 2 t 7 and 

1 

d(2) = ^ (sup|s^(<j)).- s^l + sup|s^’ (<f) |) < “ (2.2) 

i=0 

where the supremum is taken over all points of 0 < (j) < 2 t 7 where s^ *(<?>) 
and s^' ((j)) are continuous. 

We notice, in particular, that conditions a), b) and c) ensure that 

71, Tig and 71^ are square integrable on 2. We also require that the function 

w for the container walls (1.4) is sufficiently smooth in some neighborhoods 

of s = s and s = s. . We will further denote the closed rectangle s < s < s_ , 
o 1 ^ o - - 1 

0 < (|) < 2t 7 with E and use the notation 
— - o 

^77^^(<f») 

(f>g) “II f (s,^)g(s,^.) dsdcj) (2.3) 

Jo Js^((j>) 

llfll^ = 


(f,f) 


( 2 . 4 ) 



f (s,4>)g(s,(J)) dsd(j) 


(2.5) 


(f,g) 


o 



I! £11 


2 

o 




( 2 . 6 ) 


for real-valued functions f and g for which the integrals exist. 

The following lemma will be useful in the y 0 case: 

Lemma 2.1 : Let vi(s,<})) satisfy the requirements a) - c) above on the rec- 

tangle Then there exists a constant 0, depending only on s^ and s^ 
but not on ]i, such that for all r > 1, 

{li(s + ii(s ,(j))^ld^ < 3(r“^|lu 11^ + rl!y!l^) (2.7) 

iQ SO o 



Proof : For any continuous function y(x) with a square integrable first 

derivative in a < x < b and any r > 1 it holds 


|y(x) 


2(b-a) 
- 2r+l 


j: 


y'CO'^dt + 


2r 

(b-a) (2r-l) 


f. 


y(t)^dt 


in a ^ X < b . (2.8) 

((2.8) is the Sobolev inequality (see e.g. Agmon [6], p. 32) in a special case. 
The constants appearing on the right may be derived by elementary methods 
(Kato [1], p. 192-193)). Hence, we obtain for any y(s,(j>) with the assumed 
properties, any r > 1 and any 0 < <{) < 27T, 

2 7 2(s -s ) r®i ^ 2 

y(s ,(f)) + y(s,,(j))^ < y_(s,«|)) ds + — I ds 

° ^ ^ Js ^ 1 oJs 

o o 

The lemma follows by integrating this inequality over 0 < <|) < 2 tt 


Now denote the quadratic form appearing in the first Inequality in 
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Cor 2.1; Let y(s,(J>) satisfy the requirements a) - c) in the rectangle E 


and assume that, in Cl*l^)> A > A(s) > A . > 0, B > B(s) > B . >0, 

’ max - - mxn ’ max - - mxn 

C > (c(s)l holds in s < s < s, . Let Q (y) be the quadratic form defined 
max o _ _ X ^ 


above. 


Then there exist positive constants K^, K^, L^, L^, depending on s^. 


S-,, A . , 
1 mxn 


max 


B 


mxn' 


B 


max 


and C but not on y, such that 

-rr» V ‘ ■* 


max 




(2.9) 


Proof : The statement follows by using the mean value theorem and (2.7) (with 

a sufficiently large r to obtain the right inequality in (2.9)). 


We remark, that the A(s) , B(s) and C(s) in (1.14) will satisfy the above 
assumptions in some closed interval containing in its interior, because 

of (1.11) and the assumption that the unperturbed surface does not intersect 
the z-axis. 

Our sufficiency proofs will rely mainly on the following perturbation- 
result for the quadratic form Q (y) : 


Lemma 2.2 : Let A(s), B(s) and C(s) satisfy the requirements in Cor. 2.1, and 

let 6A(s,(|),), 6B(s,({)) and 6c(s,^) be boimded and integrable in E^. Further, 
let 6a^((|)) and 6a^((j>) be bounded and integrable in 0 < 4> 27T, let e > 0 and let 

^ be the functional 


# = $(6A, 5b, 6c, daj^, e) 


inf I f f |(A + 6A)y^ + (B + 6B)y^ + (C + 6c)y^| dsd<j) 

^•^0 •'s 

o 

J |(a^ + 6a^)y (s^,(j))^ + (ux 6a2^y d^| 


( 2 . 10 ) 


over all y satisfying conditions a) - c) above and such that 
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Hull = 1 <2.11) 

o 

■■ 

. |(A,V)^| < efllpjl^ + . (2.12) 

Let Ilf II denote, as usual, the supremum of jf| over the domain of f and 

= *(C.0.0.0.0.0, > 0. Ihaa poa.l« conalan. 

K and 6 independent of fiA, 6B, 6C, 6a ,6&^, and e, such that 
o O J. 


f > 4) - K 
- o 


(2.13) 


holds true, if only 

6^ = ll6Ali^ + ||6B|I^ + ||6C||^' + I|6a^|[^ + |[6aJ|^ + e < 6^ . (2.14) 

Remark ; K and 6 may depend on A, B, C, a^, a^, and s^. A similar upper 

hound also exists, and the results apply, with obvious modifications, also to 

the case ^ <0. 
o 

Proof of Lemma 2.2 ; Let Q (ji) be defined as in Cor. 2.1 and denote, similarly, 
the quadratic form in (2.10) with (Q^ + 6Q^) (}i) . Then, by (2.7) and (2.9), 
it holds for any y satisfying (2; 11) and conditions a) - c) above, 

(1 - • <S^)Q^(y) - * 5^ < (Q„ + 5Qo)(]i) < <1 + ”l ’ -'^A 

(2.15) 

where M^, M^, and are positive constants, independent of 6A, 6B, 6C, 

6a^, 6a^^, and e. Since by (2.12) f is non-increasing as function of e, 
we obtain from (2.15) 

$ < (1 + ^1 ’ ^A * (2.16) 

Now let e.g. ~ assume that y is any function satisfying 

a) - c) , (2.11) and (2,12) such that 
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(Q„ + 6Q^(u) < (1 + ^ (2.17) 

o o 

Then by (2.15) and (2.9), 

"^llo + "VoiVo + *'2 

for some positive constants ^2’ of 6A, 6B, 6C, 6a , 

6a^ and e. Obviously, then, when forming the infimimi in (2.10), only the 
subset of ]i satisfying (2.18) need to be considered. Now put 

f*(s,^) = (s - s^)(s^ - s) (2.19) 

U = G^(U - G^f*) (2.20) 

where and are chosen so as to make y satisfy (2.11) and (2.12) with 
e = 0, i.e. : 


„ 1 

^2 ~ (A,f*) » *^1 lly - C.f*ir ' 

o z o 


( 2 . 21 ) 


Notice that (A,f*)^ is a positive constant, dependent only on A, and 
s^. By (2.11), (2.12), (2.14), (2.18) and (2.21), 




i 1 - Vi 


C2.22) 


Here, and below, C^, k = 3,4,5... will be some positive constants, independent 

of 6A, 6 b, 6c, 6a^^ and e. Now form,noting that f*(s,(j>) vanishes on 

s = s and s = s- , 
o 1 ’ 


- C2i*) - Q^(li) 

n "l 

|A(s)y^fJ + C(s)yf*}dsd(J) . (2.23) 

U 
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Using the mean value theorem and Schwartz inequality in the last term, we 
obtain by (2.22) and (2.18) 


•Q^CVL - ) - Q„(M:) < 056^ + -CgiS^ 


(-2.24) 


We can now form 


- (Qo - ^o® ■ ^*^0 

< cJq^(U “ C2f*) - (1 + M^6^)Q^(y) + (2.25) 

• cJ}Q„(li - Cjt*) - Q„(!i)( -KcJ - 1 - , 


where we have used, in turn, the definition of and (2.15). Then we use 
(2.22), (2.24), (2.15) and (2.17) to obtain 


*■» i ®o + (li) + • «A 


(2.26) 


it only < 6^, where e.g. 6^ = 
independent of 6A, 6B, 6C, 5ot^» ^ 

follows immediately from (2.26). 


and K is some positive constant. 
The statement of Lemma 2.2 


For the case y® 0 we will need a slightly modified version of the 
above result; 


Lemma- 2.3 ; Let the assiamptions of Lemma 2.2 hold and let $ = $(6A, dB, 6C) 
be defined as in (2.10), (2.11), (2.12) but under the additional condition; 

U(s^,«j>) = p(s^,<j)) = 0 , 0 < 4> < 2 tt (2.27) 

Put = #(0,0,0). Then (2.13) holds true, where now 


(2.28) 
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Remark: The constants K, 6 may possibly be different from those applying to 

o 

the case of Lemma 2.2. 

The proof of Lemma 2.3 is similar to that of Lemma 2.2 and is omitted. 

We notice, though, that when the boundary terms are lacking in the qxiadratic 
form and ^ simple consequence of the mean value 

theorem and the strict positivity of A(s) and B(s), so that Lemma 2.7 is not 
needed for the proof. 
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3. SUFFICIENT CONDITIONS FOR THE CASE y > 0 
The following theorem specifies sufficient conditions for stability in 
the general case Y > 0: 

Theorem 3.1 ; Let y > 0 in (1.1) and assume that the perturbed surface 

satisfies the Euler equations (1.11) with the associated boundary conditions. 

Let further > 0 in (1.18) and let d(n) be defined as in (2.1). Then 

there exists a constant d >0, such that in (1.5) 

o’ 

6e(h) > 0 with equality iff t] = 0 in 2 (3.1) 

holds for all u satisfying the voliome constraint 6 V(ti) “ 0 in (1.6), the 
conditions a) - d) of Section 2 and 

d(n) < d (3.2) 

— o 

Proof ; 'Denote again the quadratic form in (1.14) with Q^(il)- Then by the 
definitions (1.15) - (1.19) and by the representation theorem for quadratic 
forms in Hilbert space (Kato [1], p. 322-323): 

inf Q^(y) = > 0 (3.3) 

o — o 

where the infimum is taken over all y satisfying a) - d) of Section 2 and 

llyll^ = 1 (3.4) 

(A,y)^ = 0 . (3.5) 

Now consider the expressions (1.5) and (1.6) for 6 E(ti) and 6V(ri). Firstly, 
since R(s) is positive in some open interval containing [s^,s^] there exists 
by (2.2) some d^ > 0 such that if d(2) < d^, then R(s) > R^ > 0 in E. By 
(1.8) then, when d(S) < d^ there exists some constant d 2 > 0, such that the 
functions f.CjtljS), f (n,s) and f (ri,s) are analytic in the arguments y, y 

A g V S 

and y^ in the region jy[ fly^l + [y^j < ^2 at all points of E. By (2.1) then 
at all points of E where y^ and y^ are continuous, we obtain by Taylor 
expanding at = ^: 
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9f 

fgCjljS) = g^(^,s) 


^yCn.*®) (^jS) 


n + Y + \(n>s) 

rt + I 

|a%^(0.®>n + h^.(Tl,s) 


<3. 6) 


where, for some positive 


uuii.0 Udii I- 


g 


[hA(n,s)| < Ca • d(ri)(T1g + + T]^) 

l^g(a»s) I 1 ^g ■' d(ri)n^ 


i 1 S * 


(3.7) 


if only d(Ti) < Secondly, since the perturbed surface and the container 

wall intersect at s = s.(<J)), we have w(s.(<J>)) = r)(s. (^) ,(j)) , 0 ^ <j> < 271, 

^ ^ 1+1 ^ 

i=0,l. Hence, since by (1.11) w’ (s^) = (-1) tan Y 0, we obtain by the 
inverse function theorem 


n (s •((()) ,(}>)■ 

- Sj. = -- -77^ ) ~+ \(Tl(s^ (<}>),(}))) (3.8) 

where, for some positive constants d~ and C , [h (T])j< C r|^ when jr)[ ^ d». 
Now inserting the expression (3.6) for into (1.6) and using the volume 
constraint 6V(n) = 0 we get 

3f„ 

■^^(a,s)Tldsd0 

J<^) 

.27T Sj^((i)) 

i "ii fy (o.s)a + iiy(a>s)} dsd<j) 

0 *^s^((f)) 



2it 

-iJ, 


f^(w,s)dsd(l), 


0 JA<^ 


(3.9) 
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Similarly inserting (3.6) for f. and f into (1.5) and using (1.9) and (3.9) 

A g 

gives 

2ir s^((j)) 

6E(n) := I I + fg^^(0,s) - Xf^^^(£,s)}n + h(Ti,s)|dsd(j) 

■ •'o •'s (4>) ~ ~ 

O 

|cosYf^(w,s) - f^(^,s) + fg(w,s) - Xf^(w,s)Jdsd^ , (3.10) 


where h(_g,s) = h^(r],s) + h (ji,s) - Xh (jj,s) satisfies by (3.7), 


lh(n,s)j < C ' d(n)(rij + 0^ + 0^) 


(3.11) 


at all points of Z where 0^ and are continuous. Here again C is some 
positive constant independent of p, and we have assumed d(ll) < as in 
(3.7). Introduce for the second integrand in (3.10) for convenience 


fj^(s) = cos Y f^(w,s) - f^(^,s) + f<g(w,s) - Xf^(w,s) 


(3.12) 


Then, by (1.7) - (1.9), fi,(s ) = f (s ) - 0. Assuming w(s) to be sufficiently 

it O K J. 

smooth for f to be twice continuously differentiable in some neighborhoods 
K. 

of s = s^ and s = s^ and using the definition (1.7) for A(j), Taylor's theorem 
gives 


I 


fj^(s)ds = 


i+1 


1=0 


V(Si)(s.(<^)-s.)^ 


+ hj^((j)) 


(3% 13) 


where, for some positive constants d^* and independent of <j), 
ih^(<{>) I < C^|s^(<{)) - if only d(Z) < d^. Using (3.8), (3.13) and 

the notation in (1.13) and (1.14), (3.10) takes the form 
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2w s^(<J>) 

6E(n) =j! J I + C(s)n^ }■ + h(j^,s)|dsd(J) 


2it 


+ Sa^((f))} n(s^(<{>) ,4>)^ + {otj^ + (<|))}ti(s^ .(<)>) >^)' 


d(J) 


(3.14) 


where (3.11) and 


6a^(<!)) I < Cj^ ' d(n) , 1=0,1 


(3.15) 


hold, if only d(T]) + d(S) < d^, where again and d^ are positive constants 
independent of n. Treating the last integral in (3.9) by Taylor expanding 
the integrand in the same way as in (3.12) - (3.13) and- using (3.8) the volme 
constraint gives 

27T 

I ( I A(s)ndsd^[ 

•'O *^s 


2tt s^ 2ir 

- ^if f Tl(s,(j))^dsd(f) + K^2 r {n(s^((j)),(j))^ + n(Sj^(cf))„(f))^} dct) 

•^0 *^s •'O 

0 


(3.16) 


if only d(n) + d(S) < d^, for some positive Ti-independent constants 

and d_. 
b 


In order to relate the above expressions for 5E(n) and 6V(n) to the condi- 
tions (3.3) - (3.5) (which involve integrals over the rectangle S ) we now perform 
a transformation of variable in (3.14) and (3.16); 


S ^ " s 

s = s(s',(})) = s (({)) +- — (s (<j)) - s (<j))) (3.17) 

O S- S J. o 

1 o 


iJ(s.’,<f>) = n(s(s’,^),4>) 


(3.18) 
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We note, in particular, that p will satisfy the. requirements a) - c) of 
Section 2 if only dCS) is sufficiently small to ensure I s, (<{») - s (<|)) I > s . > 0 

1 O' - mm 

in 0 < ^ < 2ir. Furthermore 


s = s’ + k^(s’,<})) 




r)g(s(s’ = (1 + 


(3.19) 


ri^(s(s‘ ,(|)) ,(|)) = p^(s’,cf)) +k^(s’,4>) • p^,(s’,(|)) 




where, for some positive constants and independent of T|, 


sup •{ lk^(s' ,(f)) [ , |k 2 (<{>) Mk^Cs’ ,4>) 1} 5 • d(S) (3.20) 


at all s^ < s’ < and all 0 < (J) < 2 it where s^((})) and s^(4>) are continuous, if 
only d(2) < dg. 

Inserting (3.17) - (3.19) into (3.14) and (3.16) we get: 

2ir s^ 

'SE(ti) = I I I {(A + 6A)Pg, + (B + 6B)p^ + (C + 5C)y^}ds’d(}) 

*^0 *^s 

o 

2ir 

+ ) ^%(<i>))vt(s^,({))^ + (a^ + 5a^(c|)))y(s^,<|))^}d4> (3.21) 

•'o 


and 



(A + 6’A)pds'd$j 


2ir 

< K^3llpfl^ ^ ^2 J + U(s^,(j))^fd<j) , 


(3.22) 
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where, by (3.11), (3.19) and (3.20), 6 a, 6b, 6C and 6'A are functions 
of s' and ({>, integrable in and such that 

sup{| 6 A(s',(i))| + | 6 B(s‘,<|))j + l6C(s',4>)| + j 6 ’A(s’,(|))|} < C 2 {d(n) + d(E)} (3.23) 

when d(n) + d(Z) < d^ where C 2 and d^ are some positive constants, independent 
of r|. The supremum in (3.23) is formed over all points of 2^ where the functions 
are continuous. Note, in particular, that the argument in the functions A, 

B and C in (3.21) and (3.22) is s', whence in obtaining (3.23) we have made 
use of the fact that A, B and C of (3.14) (as defined in (1.14)) have a 
bounded derivative in 2 (they are in fact analytic ftmctions of s),, if only d(2) 
is sufficiently small. 

Using Lemma 2.1, (3.23) and the mean-value theorem in (3.22) gives in a 
straightforward way 

|(A,u)^| < C 3 • {d(Ti) + d(2)} (llyjl^ + + lliill^) (3.24) 

if only d(Ti) + d(2) < d^, where C_ and d„ are some positive constants, 

— o o o 

independent of r|. 

Now using Lemma 2.2 on (3.21) and (3.24) (noting (3.15) and (3.23)) we 
conclude by (3.3) that, e.g., 


6E(n) > (3.25) 

if only 

d(u) + d(2) < d^ (3.26) 

for some ri“ii^dependent constant dg. Since, however w(s^($)) = Ti(s^((|)) ,^) , 

0 < (j) < 2ir, i=0,l, it holds 


nrf,(s. (<{)), <j>) 

s'fd)') - X 3r. ■ 

i^^ • w'(s ((^)) - ri„(s ((})), 4») 

X 5 X 


(3.27) 


at all points of ■ the boundary of 2, 


where U 


<!> 


and T) 

s 


are continuous. 


Since 
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x4*l 

further w’ (s^) = (-1) tan Y ^ 0, it follows from (3.27) that 

d(S) < C^d(ri) (3.28) 

if only d(ri) < d^^, where d^^ and are some positive, ri-independent constants. 
Hence, by (3.26) and (3.28), (3.25) holds true, if only d(T]) < d^, where 
d^ is some n- independent constant, and the first part of the statement (3.1) 
of Theorem 3.1 follows. Since, by (3.18), IIt)1I = 0 if and only if 
the second part of the statement also follows from (3.25). 
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4. SUFFICIENT CONDITIONS IN THE' CASE Y =‘ 0 
In this section we will show that the statement of Theorem 3«1 also holds 
in the limit Y = 0 (defining the limit of the boundary, conditions in (1.17) as 
|a. 3 = 0,1 to be (1.20)), if the curvatures of the imperturbed surface and 
the container wall' are not equal at the unperturbed- contact lines . The 
statement is, in more precise terms, the following: 

Theorem 4.1 ; Let Y = 0 in (1.1) and assume that the unperturbed surface 

satisfies the Euler equations (1.11) and the associated boundary conditions 

at s = s and s = S- Assume further that the function w(s)- of (1.4) 
o 1 

satisfies 


w"(s^) <0 , i=0,l (4.1) 

Let X be defined as in (1.17 - 1.19), but with the boundary conditions 
o 

of (1.17) replaced by (1.20), and assume that > 0. Let d(ri) and d(2) 

be defined as- in (2.1) and (2.2). Then there exists a constant d >0, such 

o 

that in (1.5) 


6E(ri) > 0 with equality iff ri = 0 in Z (4.2) 

holds for all D satisfying the volume constraint 6 V(ti) = 0 in (1.6), the 
conditions a) - d) of Section 2 and the condition 

d(n) + d(E) < d^ . (4.3) 

Remark ; In the case Y “ 0, both the energy (1.1) and the liquid volume of any 
configuration are unchanged if the liquid-vapor -interface is continued by "wetting" 
the "dry" container walls (and we will make use of this property of the "wetting" 
perturbation in the proof below). With E in (4.2) (see also the beginning of 
Section 2) for such "wetting" perturbations is meant a closed domain, con- 
taining no interior points where the liquid-vapor interface coincides with the 
container walls (i.e. possibly wetted portions of the container walls are 
not included in the E). Note, further, that since by' (1.11) w’ (s^) = 0, 1=0,1, 
and since w(s) < 0 for s^ < s < s^, we must have w"(s^) < 0, 1=0,1. w" (s^) = 0 
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means by (1.4) that the curvatures of the unperturbed surface and the 
wall coincide at the unperturbed contact' line s = ( and is a limiting case 

among all solutions to ( 1 . 11 ) that correspond to permissible liquid-tank 
configurations) . 

Proof of Theorem 4.1 ; Let again Q (jl) be the quadratic form, in (1.14). 

Then, as in (3.3) 


inf Q (y) = X > 0 (4.4) 

o — o 

where the infimum is taken over all y satisfying a) - d) of Section 2 and 
the conditions: 


llyll = 1 
o 


(4.5) 


(A,y)^ = 0 


y(s^,<{)) = y(s^,^) = 0 , 0 < (j) < 2 tt 


(4.6) 

(4.7) 


Here we denote, for brevity. 


f 5 .(H,s) = f^(Tl,s) - f^(0,s) + fg(jl,s) - Xf^(n,s) 


(4.8) 


where f^, f^ and f^ are the functions in (1.8). Then, as in (3.10), we have 
since Y = 0 , 


■cr 


6E(ti) = 


fg(] 1 ,s)dsd(|) 


0 *'s^((f>) 


-n 

0 


fg(w,s)dsd<{) 


(4.9) 


Now denote, for convenience in the following. 


Z 4 . = complement of (SOS^) w.r.t. Z 


Z_ = complement of (E<^Z^) w.r.t. 


(4.10) 
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Let ri(s,^) be the function obtained by extending ri(s,(|)) by "wetting" 
those parts of the container walls which, were "dried" because of the pertur- 
bation n, i. e. : 


ri(s,(j)) = 




Tl(s,^) in Z 


[w(s) in Z- 




Put, further,' ZuZ = ZuZ_ = 


Z' uZ , with the boundaries 
o + 


s^((}») = min|s^((f)) 


s^((j)) = max|s^(c|>) 


0 < (j) < 2ir . 


(4.11) 


(4.12) 


Then it follows from (4.9) (and, more generally, from the remark after the 
statement of Theorem 4.1), that 


2^ Ji(<j)) 


5E(n) = 


u. 


'0 J s ((j)) 
■ o • 


fE(5..s)dsd<|) 



fg(w,s)dsd<j) 


(4.13) 


where A<{) is defined from s (<j)) and i- (ij)) as in the last formula in (1.7) (in 

o ^ ^ 

other words, the domains of integration in (4.13) are Z and Z+, respectively). 
Put further. 


0 in Z 


h*(s,(|.) = 


\w(s) in "Z..J, 


V(s,{j>) = ri(s,(f>) - •n*(s,<j)) in Z 


(4.14) 


(4,15) 


Then, noting that f^(^,s) = 0, it follows from (4.13), that 
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2it 

6E(Ti)=j I |fg(n* + v,s) - f^(Ti*,s)(dsd(i) (4.16) 

*'0 


-Similarly-, we obtain by (-1.6) 

f 2ir s^(<i>) 

I ]f^(Tl* + V,s) - f^(jn*,s)}dsd(|) (4.17) 

0 *^s^((})) 


We note, in particular, that ri, r)*, V and s., i=0,l will satisfy the require- 

^ 2 3 

ments a) - d) of Section 2 and further that, since- w(s) = w"(s^)(s-s^) /2 + 0(s-s.) 
in the neighbourhood of s = s^, i = 0,1, it holds 


d(S) < d(E) 


d(n*) < c^d(E) 

d(fi) < c^]d(n) + d(E)} 
d(v) < c^|d(n) + d(E)[ 


(4.18) 


if only d(ri) + d(S) < d^', where and d^ are some rj-independent constants. 
Also, as in (3.6) the integrands in (4.16) and (4.17) are analytic in the 
arguments q*, ri*, ri| and V, v^, in some region [t]*| + |^[ < where d 2 
is an Ti-independenf constant. By (4.18) then, we may Taylor expand the 
integrands in (4.16) and (4.17) at all points of E where 71* and-^ are contin- 
uous, if only d(n) + d(E) < d^, where d^ is some Tl-independent constant. 

We obtain 


” fg(jl*,s) 

= WH*»s)^ • V+ |v'^fg^(0,s)V + hg(T]^*,V,s) 

and 

3f„ . 

fy(ir + ^,s) - fy(ll*,s) = ^(0,s) • V + h^(Tl’^,V,s) 


(4.19) 


- C4-.20) 
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Where, by (4.18) and noting from (1.9) that fg^(£,s) =.0^, 

- I . 

]h^(Ti*,V.s)l < C 2 |d(Ti) + d(Z)[(vJ + V^+ v2) (4.21) 

|h^(Ti*,v,s)[ < C2|d(Ti) + d(Z)} (|Vg| + |v^| 1- |v|) (4.22) 

if only d(ri) + d(S) < where d^ and are some D-independent positive 
constants. 

Now, by (1.8), (1.9) and (4.8), the first term to the right in 
(4.19) is of the form 




0 in E 

o 


= |a(s) + D^(s)}w* (s)Vg + |c(s) + D^(s)}.w(s)V in S+ (4.23) 

where A(s) = OR.(s) and C(s) are the functions introduced in (1.14) and where 

Cs)^}- (4.24) 

if only d(E) < d^, where d^and are some positive, ^-independent constants. 
Further, I^^(s) is continuously differentiable in the interior of E^ if w"(s) 
is continuous there. 

Now consider the contribution of the terms in (4.23)' to the energy 
integral (4.16): 



(4.25) 
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By partial integration of the first term over A(j) = (s^(<j)) ,s^)u(s^,s^ ((}>)) , 
noting that w'(s.) = 0, i=0,l, and that V(s .(<}>) = 0» 0 < <f) < 2ir by (4.15) 

X X * 

we get 


Jo Js Of') 

O 

2ir 


fg^( 2 i*,s)^ • V dsd(|) 


(4,26) 


n |-|^{(A(s) + D^(s))w'(s)} + |c(s) + D^(s)}w(s)|vdsd(f) 

A 


'0 *'A(j) 


Note that it holds in A<f), if d(2) is sufficiently small. 


A(s) = OR(s) > R > 0 
— -o 

w(s) = w"(s^)(s - s^)^ + 0(s - s^)^ , 

v(s,(f)) > 0 (with strict inequality - in interior points) , 

(4.27) 

i (d>) < s and s- < s, (d)) , ' 0 < A < 2ir . 

0-0 1-1 -“ 


The last two formulas follow from (4.15) and (4.12), respectively. Using 
(4.27) and (4.24) we then see that since w"(Sj^) < 0, i=0,l, the last 
integrand in (4.26) is non-negative in Aip (with -A(s) • w"(s) • V(s,4>) as 
the dominating term), if only d(2) is sufficiently small. Consequently, 


2ir 

U 




,s) * V dsdcj) > 0 


0 •'Acf) 


(4.28) 


if only d(E), < dg where dg is some u- independent constant. (4.16) and (4.19) 
then give, using the notation of (1.14) 


21T i^((f)) 

6E(n) > f I I {A(s)Vg + B(s)v^ + c(s)v^ + 2hg(Tl*,V,s) } dsd<t> (4.29) 

•^0 *^%(<!>) 
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where V(s,i) satisfies the volume constraint 6V = 0, which by (4.17) and 
(4,20)' is of the form 


hy (jL* s ) ds d({> 


Sl((f>) 

2tt 1^(4)) 

1 A(s) • vdsd^ = -cr • 

f [ 

•^i^((j)) 

“'0 *^s^((f)) 


(4.30) 


Again, the functions h^ and h^ of (4.29) and (4.30) will satisfy (4.21) 

and (4.22), if only d(n) + d(S) < d^, where d^ is some ^-independent constant. 

In order to relate (4.29) - (4 ..30) to the properties (4.4) - (4.7) of the 

atic form Q , i 
o 

change of variable: 

s = s(s’,^) = s ('(j)) + ^(s ((J>) - s ((j))) 

O S- “■ s J. o 


quadratic form Q^, we proceed as in the proof .of Theorem 3.1 by introducing a 


s - s 


(4.31) 


ly(s’,4») = v(s(s',(J>),(J)) 


(4.32) 


(4.31) takes t onto the rectangle 2 and the transformation satisfies (3.19) 
(with n replaced by V) and (3.20). Furthermore, by (4.15) and (4.32), 

U(s ',<{)) will satisfy (4.7). Using, in addition, (4.21) and (4.22) we obtain 
in the same way as for (3.21) and (3.24), 


27T s 


5E(n) = f j* |{(A + 6A)yJ, + (B + 6B)ii^ + (C + 5C)u^}dsd^ 

•'O *'s 


(4.32) 


where 6A(s',<j>), 5b(s’,4)), Sc(s',(f)) are bounded and integrable on with 


|(A,y)^| < C^jd(n) + d(Z)}-(l|yg,||^ + Ily^li^ + llyll^}, 
!I6 aI|.^ + I16BII^ + llscll^ < c^{d(ri) + d(Z)} 


(4.33) 

(4.34) 


if d(ri) + d(Z) 5 *^8’ some positive constants, independent 

of U. We can then apply Lemma 2.3 to conclude by (4.4), that e.g. 


6E(ri) > -^lyllj 


(4.35) 
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if only d(Ti) +’d(S) < d^, where d^ is some r)-independent, positive constant 
Since, furthermore, by (4.11), (4.14) and (4.15) V = 0 in S if and only if 
r| = 0 in E and since by (4.32) Upll^ = 0 if and only if v = 0 in Z, (4.35) 
completes the proof of Theorem 4.1. 
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5. EXPLICIT FORM OF THE JACOBI EQUATIONS. CONCLUDING REMARKS. 

The functions A(s) , B(s) and C(s) of (1.14) are found by straightforward 
manipulations from (1.8) to be (we put R = R.(s), Z = Z(s) for brevity) s 

A(s) = CTR(s) 


B(s) = cr/R(s) 


C(s) = -20R" + pg]RR' - ZZ’ + ZR(R"Z' - Z’'R* )[ 


+ -X]Z’ - R(R"Z’ - Z"R’){ 


(5.1) 


Using (5.1) the Jacobi differential equations (1.17) become (putting 
B^ = pg/ff and = A/a as in (1.11)): 

RU" + R’UL - - 2R" + B ^RR’ - ZZ’ 

xk ik (R o' 

+ ZR(R"Z' - Z"R’)} + H^|Z' - R(R"Z' - Z"R’ )[| 

" ^ik^ik ; k= 0,1,2... 

i = 1,2,3... 

s < s < s. , (5.2) 

o — — 1 

with the botindary conditions 

(-l)'^tan y y’(s^) = -{sin^Y(R"Z' - Z"R’) 

i=0,l. (5.3) 


- cos^yw"[g^g . y(s^) 
i 


Note that (5.3) tends formally to p(s,) - 0 as y ^ 0, if w"(s.) ^ 0, i-0,1 , 

X 

which is the case covered by Theorem 4.1. 

Finally we remark that the conditions b) and d) of Section 2 may probably 
be slightly relaxed by introducing concepts from the theory of Lebesque 
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integrals. .With regard to the physical background of the problem, however, 
such extensions of the set of permissible perturbations do not seem very 
meaningful. 
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STABILITY OF AXISYMMETRIC, ANNULAR FLUID INTERFACES 
AT ZERO CONTACT ANGLE 


Ilkka Karasalo 


ABSTRACT 

We study the stability, in terras of minimal total potential energy, 
of liquid configurations in axisymmetric containers, such that the 
liquid-vapor interface is annular and meets the container walls at zero 
contact angle. The proper limits of sufficient and necessary conditions 
for stability, respectively, as the contact angle tends to zero, are 
formulated in terms of the Jacobi accessory differential equations. The 
stability is shown to depend crucially on whether the equilibrium liquid- 
vapor Interface stays inside the container or not when continued 
analytically past the three-phase contact lines. 
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1. INTRODUCTION 

We shall study in this paper the stability of certain configura- 
tions of liquid partially filling an axially symmetric tank in a 
..gravitational field directed along the axis of symmetry. We require, 
that the tank shape and the liquid volume are such that the liquid- 
vapor interface is annular, i.e. it does not intersect the axis of 
symmetry, cf. figure 1; 


Container 

wall 


XBL 766-3054 

Figure 1 : Example of permissible liquid-tank configuration 

and associated coordinate system. 
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A configuration is in stable equilibrium if and only if it strictly 
minimizes the total static potential energy of the system. 


E = a • (A, - cos Y ■ A ) + E 
r w g 


( 1 . 1 ) 


among all nearby configurations with the same liquid volume V. Here 

a > 0 (the liquid-vapor surface tension) and 0 < y < tt (the contact 

angle between the liquid-vapor surface and the container wall) are 

constants, A^ and A^ are the areas of the liquid-vapor and the liquid- 

wall interfaces, respectively, and E is the gravitational potential 

S 

energy of the liquid. This constrained minimization problem has 
received much attention in the literature, see e.g. Huh [5] and 
Gillette [4] for extensive lists of references. By a suitable choice 
of variables, it may be viewed as a variable-endpoint problem of 
variational calculus .([4] p. 21 and p. 145). When Y > 0, this approach 
results in conditions which distinguish between stable and unstable 
cases in a rather satisfactory way. There appear to- be fewer rigorous 
results, however, concerning to what extent these stability conditions 
also apply to the limiting case Y = 0 (cf. [4], p. 23). The purpose 
of this paper is to analyze this limiting case for axially symmetric 
liquid configurations of the above kind. More specifically, we shall 
look at necessary and sufficient conditions, respectively, for 
minimum of E based on the Jacobi accessory minimization problem for 
the second variation of E (see e.g. Akhiezer [1], p. 69), as Y 0. 

The formal limits, as y Oj of the boundary conditions associated 
with the Jacobi accessory differential equations depend crucially on 
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whether the curvatures of the equilibrium liquid-vapor Interface and 
the cp.nt.ainer wall, coincide pr no.t at the three-phase contact lines. 

In the latter case these limiting" boundary conditions will be of the 
fixed end-point type (when using a parametric representation of the 
surfaces, see further , Section 2.1 below). 

We will show, firstly (Theorems 3.1 and 3-2) that the stability 
conditions (sufficient and necessary, respectively) based on the fixed 
end-point boundary conditions in fact apply to (1.1) with y = 0 if 
only the analytic continuation of the equilibrium liquid-vapor interface 
does not penetrate the container walls at the three-phase contact lines . 
Secondly (Theorem 3.3), we show that if the analytic continuation of 
the equilibrium liquid-vapor interface does penetrate the wall at either 
of the contact lines, the configuration will be unstable regardless of 
the conditions on the second variation of E. 
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2. NOTATION AND SOME PRELIMINARY RESULTS 
2.1 The Euler-Lagrange and Jacobi Conditions 

It will suffice to consider (1.1) at small perturbations from 
axially symmetric configurations. We will use a parametric arc-length, 
normal displacement representation of the surfaces (see e.g. Reynolds, 
Saad, Satterlee [8]), Thus the unperturbed liquid-vapor interface is 
described by 


r = R(s) s < s < s, 

o - - 1 

z = Z(s) 0 < (|) < 2 tt 


( 2 . 1 ) 


in the polar co-ordinate system of figure 1, where s is the arc-length 
along the curve of intersection between the interface and any plane 
4) = constant. Then the equations 


r = R(s) - n(s,(i))Z' (s) 
z = Z(s) + ri(s,({))R' (s) ‘ 


s (4>) < s < s (0) 

O — — 

0 < <J) < 2ir 


( 2 . 2 ) 


describe a surface obtained by moving each point of the surface (2.1) 

the distance r|(s,(i)) in the direction of the normal at (s,4>). (In 

general, since we want the perturbed surface (2.2)' to intersect the 

container walls, the functions R and Z of (2.1) must be continued to 

some open interval containing [s^,s^]. A convenient way of doing this, 

which we will use in the sequel, is provided by the differential 

equations (2.9) below). Similarly, in some neighborhood of the 

unperturbed contact lines (s = s and s = s in (2.1)), the container 

o 1 
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wall will be given by 


r = R(s) - w(s)Z'(s) 

; 0 < (f) < 2ir. (2.3) 

z = Z(s) + w(s)R’ (s) 


Then, denoting by 6E(n) and 6 V(ti) the increments of the energy (1.1) 
and the liquid volxime at the perturbation (2.2), we obtain in a 
straightforward way 


2tt s-(<j>) 


6e(ti) 


✓ J. V 


0 


2t{ 

- I 


cos Y ' f.(w,s) - dsd<^ 


0 A({) 


2ir _s, (({)) 


6V(n) 


‘'o *'s^(4>) 


f^(_n,s)dsd^ 


-// 

•'rt A J 


fy(w,s)dsd<{) 


0 


(2.4) 


(2.5) 


where we have put 

ri = Ti(s,(f)) = (n(s,(}!),Tig(s,(|)),n^(s,4)))^ 

w = w(s) = (w(s) ,w* (s) ,0)^ (2.6) 

A^ = the interval (s^((|)),s^)u(s^,s^(<j))) 


and, denoting R = R(s), Z = Z(s), 
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^ ~ nz')^(ng + (1 + n(R"z' - z”r'))2) 

+ n?(i + n(R"z’ - z"R*))^f^^^ 

<p 

f^(Tl,s) = n|l +-^(R"Z’ - Z"R’)f(;R - |Z') 

fgdl.s) = Pg * (Z+|r’) . f^Cn.s) (2.7) 


Here p is the constant liquid density and g is the gravitation constant 
with g > 0 if the gravitation force acts towards the negative z-axis 
in figure 1. 

The condition, that all first order Ti-terms in 6E(ri) should 
vanish for all p such that 6V(p) = 0 then leads to 



(0,s) 


3f 

+ 3^(0, s) 



= 0 


ins <s<s,, with the boundary conditions 
o — — 1 

cos Y ^®i^’®i^ ~ *®i^ ^ ’ i=0,l , 


( 2 . 8 ) 


where X is a constant (the Lagrange multiplier). Putting B = pg/P, 

o 

H = X/p, using (2.6), (2.7) and the identity R'(s)^ + Z'(s)^ = 1,(2.8) 
o 

becomes the Euler-Lagrange boundary value problem 


■R" = -Z' (B Z - H - Z’/R) 
o o 

Z" = R* (B Z - H - Z'/R) 
o o 

w(s^) = 0 

w'(s^) = (-l)^"^^tan Y 



(2.9) 
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(Without loss of generality we have excluded the contact angles 
TT - rif allowed by (2..8)),. 

Assuming now that (2.9) is satisfied, the condition, that all 
second order ri-terms should give a non-negative contribution to 
6E(ti) for all rt such that 6v(n) = 0, takes the form 


27T s 


QqCIL) J* |A(s)yg + B(-s)y^ +' C(s)vi^ } dsd(j) 


0 , s 


2-n ° ' 

I |a^y(s^,({))^ + a^y(s^,(j))^| > 0 

r\ 


( 2 . 10 ) 


for all y(s,4>) such that 



R(s)y (s,4))dsd({) = 0 


( 2 . 11 ) 


Here we have put 


^Arin^-’®) "" diag|c(s) ,A(s) ,B(s)f , (2,12) 

- ^~p" • ■fj{cos Y f^(w(s),s) - f^(0,s) + f (w(s),s) ~ Xf^(w(s),s)|g^g 


= a. 
1 


i=0,l. 


(2.13) 


By (2.12) and (2.7), the A, B and C of (2.10) are 

A(s) = OR(s) 

B(s) = a/R(s) 

C(s) = -2aR" + pgjRR’ - ZZ' + ZR(R"Z’ - Z"R')^ 
+ X|Z' - R(R"Z'- - Z"R’)^ 


(2.14) 
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By (2,9), since R(s) >R. >Oins <s<s,,A, B and C will be 

— mxn o — — 1 

smooth (in fact analytic) in some .open interval containing 

and A(s) > A . >0, B(s) > B . >0 will hold there. By standard 

— min — min 

results for symmetric, semibounded quadratic forms in Hilbert space 
(see e.g. Kato [7], p. 322 and 352-353), (2.10) may be analyzed in 
terms of the eigenvalues and eigenfunctions of an associated self adjoint 
differential operator: 

Tp = -|^(A(s)y^) - |^(B(s)y^) + C(s)y 

in < s < s^, 0 5 with the boundary conditions, that y 

should be periodic in 4> with period 2 tt and 

A(s )y is ,{^) = (-l)Vy(s ,4>) , 0 < (J) < 2it 

i = 0,1 . 

T has a complete, orthogonal system of eigenfunctions of the form 

|li^^(s) =os , ju.^(a) sin 

with associated eigenvalues (ordered increasingly in the 

index i) , determined from the boundary value problems 

- |j(A(s)y^^(s)) + |k^ . B(s) + C(s)}y^^(s) = 

A(s )y' (s ) = (-l)^a y (s.) ; j=0,l ; (2.15) 

3 J J IK J o o 

X— _L y 

k=0,l, 2, . . . . 
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We notice, that all eigenfunctions but those with k = 0 satisfy the 

constraint ( 2 . 11 ) and that, since B(s) > 0 in s < s < s, , the 

o - - 1 

eigenvalues K are increasing functions of k. It then follows that 

( 2 . 10 ) with the side-condition ( 2 . 11 ) holds for all y in the class of 

continuous functions ins <s<s_,0<(i>< 27T, which are periodic 

o — — 1 — — 

in (j> with period 2 ir and have square integrable first derivatives 
(see e.g. Kato [7], p. 322-323, Cor. 2.3) if and only if 


K 

o 


^2'^20’'^lli 


(2.16) 


where, denoting (f,g)^ 
solutions to 


•'0 *'s 


f (s,^)g(s,(J))dsd(i>, 6 and are the 





^l^'^lO’^lO^o'*’ ^2 ^^20 *^20^0 


(2.17) 


with 62 = 0 if the solutions are non-unique. 

(2.15) - (2.17) are(the equivalent of)the Jacobi accessory 
boundary value problems for our constrained minimization problem. By 
(2.13), (2.6) and (2.7), the boundary conditions are 
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Vii^CSj) = 0 , j=0',l , (2.19) 

as Y 0. 

2.2 Permissible Perturbations. Two Lemmas. 

By the contact lines of the liquid-vapor interface we mean the 

two closed curves within the container wall, any open neighborhoods 

of which intersects the interiors of both the liquid and the vapor 

inside the container (cf. fig. 1). The contact lines determine 

by (2.2) a closed region s (<j)) < s < s, (d)) , 0 < (J» < 2ir in the 

o — — i — _ 

(s,(j)) -plane. We denote this region with S and require the following 
regularity properites from Z and the associated function q: 

a) Tl is continuous in Z and periodic in <j) with period 2tt. 

b) rig and are continuous in Z except possibly at finitely 

many isolated points or finitely many piecewise smooth curves 

with finite length. In particular, ll is piecewise 
^ s 

continuous as function of s for all 0 < < 27T. 

c) d(n) = supiri(s,^)| + sup|rig(s,(j)) I + sup jri^(s ,(()) | < (2.20) 

where the supremum is taken over all points of Z where 
Tig and r|^ are continuous. 

d) s^(4>), i=0,l are continuous and such that s^(4>) - s^, 1=0,1, 
change sign at most finitely many times in 0 < (j) < 2tt. 

s^ (({)), i=0,l, are continuous except possibly at finitely 
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many points in 0 < <}) < 2 it and 
i • 

d(2) - (sup|s^((j)) - s^l + sup|s^((})) j) < « (2.21) 

i=0 

where the supremum is taken over all points of 0 < ({» < 2 it 
where s^($) and s^((|>) are continuous. 

Remark : The sufficient conditions to be considered below will 

ensure the stability of the surface (2.1) with respect to all pertur- 
bations (2.2) which satisfy a) - d) above and for which d(ri) + d(I) 
is sufficiently small. Thus in terms of variational calculus (see 
e.g. Bolza [2], p. 68-70) the extremum will be "weak”. The detailed 
assumptions imder b) and d) are introduced for simplicity in what 
follows, and could be relaxed slightly by introducing more advanced 
concepts from the theory of Lebesque integrals. With regard to the 
physical background, however, nothing essential is lost by the above. 

We will denote the closed rectangle s^ < s < s^, 0 < <j) < 27T by 

Z and use, for any f(s,(J>) which is square integrable on Z , 
o o 

llfll^ = (f,f) (2.22) 

o o 

where ( , ) is defined as in (2.17)* Then the following result will 
o 

be- useful: 

Lemma 2.1 : Let satisfy the requirements a)-c) on and t he rand i t i on 
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u(s^,4>) = y(s^,4)) =0 , 0 < <p < 2n. (2.23) 

0 (li) 2»e defined as in (2.10) and assume that A > A(s) > A . > 0, 

o — max - - min 

B > B(s) > B . > 0 and fc(s)f < C hold in s < s < s. Then 

max - - min - max o ~ - 1. 

there exist positive constants K , K,, L and L,, depending on s , s,, 

o 1 o 1 ^ o 1 

A A , B , , B and C but not on u, such that 
min max min max inax 

K 0 (y) - L llpll^ < liy II^ + IlyJi^ < K-Q <y) + L.llyli^ . (2.24) 

oo— o o— So <po-io — . J o 

(2.24) follows from the mean value theorem in a straightforward 
way and the proof is omitted. (By use of a Sobolev- type inequality 
(see e.g. Kato [7], p. 193), condition (2.23) could in fact be omitted, 
and this stronger result could be used for a similar treatment of the 
case Y ^ 0» Karasalo [6])- 

Before stating our second Lemma we need some further notation. 

For any ]i(s,<{)) satisfying a) - c) on and (2.23) we put, for clarity 

Q (A,B,C,y) = Qju) (2.25) 

o — o — 

where Q (y) is defined in (2.10). Further, if 5 A(s,<|)), 5B(s,<j)) 
o — 

and 6C(s,(j)) are bounded and integrable on and £ is a positive 
constant, we put 

$ = 4>(6A,6B,6C,£) = inf Q^(A + 6A,B + 6B,C + <SB,y) (2.26) 

over all y satisfying a) - c) on 2^, (2.23) and the conditions 

ilyli^ = 1 , 


(2.27) 
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|(R,u) I < eClIyJ^ + ilyJI^ + llyll ). (2.28) 

O ”” SO o o o 

If f(s,4>) is bounded on S , we will let, as usual, !lf|| denote the 

o ^ 

supremum of |f(s,<J>)| over (s,it)) in Then we have 

Lemma 2.2 ; Let A(s) , B(s) and C(s) satisfy the requirements in 
Lemma 2.1 and let 6 a(s,^J, and Sc(s,if>} Jbe bounded and 

integrable on Let e > 0, define $ as in (2.25) - (2.28), denote 

” ^(^f0,0,0) and put 

6 ^ IldAii + II6BII + IldCll + e . (2.29) 

GO OO 00 

Then there exist positive constants C and 6 , independent of Sa, Sb, 

o 

6c and e, such that 

1$ - $qI f ^ (2.30) 


holds true, if only 8 < 6 . 

- o 

Proof : Throughout this proof, M^, N^, 6^, 1=1, 2, 3, will 

denote positive constants, independent of 6A, 6 b, 6c, £ and y. Let 
y satisfy conditions a) - c) on (2.23) and (2.27). With the notation 
of (2.25), put for brevity, 5Q^(ji) = Q^(5A,6 b, dC,jj) . will be 

well defined because of the assumptions.) Then there exist M 2 , 
and N 2 , such that 

(1 - M,6)Q^(y) - N.6 < Q^(y) + 6Q^(y) 
i o — i — o — o — 

( < (1 + M.6)Q (y) + N,5 

> , — ±0 — i 


(2.31) 
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because of Leiraa 2.1 and the mean value theorem. Noting by (2.26) 
and (2.28), that $ is a non- increasing function of e and that, by 
(2.10) and (2.25), Q^(A,B,C,ji) is a linear function of A, B and C 
we obtain from (2.31) 


$ < (1 + M,6)$ + N,6 (2.32) 

It follows, that we need only consider those y which satisfy, e.g., the 
additional condition "h 5 2M^6)|#^j + 2N^5 when 

forming the infimum in (2.26). By (2.31) and Lemma 2.1, however, for 
all such y 

llyjl^ + !!yj^ < | + n, = m_ (2.33) 

So fo-Z'o' I 3 

if only, e.g., S < 5, = . Let y be any function satisfying a) - c) 

— 1 '^^1 

on (2.23), (2.27), (2.28) and (2.33). Put 

f*(s,<|>) = (s - s^) (s^ - s) 

y « c^(y- C 2 f*) 

where and C2 are chosen so as to make y satisfy (2.27) and (2.28) 

with e = 0, i.e. C, = (R,y) /(R,f*) , c. = 1/liy - C,f*li . Noting 
that (R,f*)^ is a positive constant, dependent only on s^, and 
R(s), it follows from (2.27) - (2.29) and (2.33) that 

IC2I < M^5 ; ^ > 1 - M56 (2.34) 

Furthermore, since f* vanishes on.s = and 5 = s^, y satisfies (2.23). 
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Then form 


= (cj - l)Q^(u) + C^C^Q^Cf*) - 6Q (y) 

X o — XXO — o — 

27t 

- I |A(s)ygf* + C(s)yf*f dsd(j> . 

“'O *s 


Here we use Lemma 2.1, (2.31), (2.33), (2.34), the mean value theorem 
and the Schwartz inequality to find upper bounds for the terms to the 
right. We obtain, that for some Mg, 6^ 




holds true, if only ^ < <52* Hence 


$ _ $ < M, ♦ 6 

o - 6 


(2.35) 


(2.36) 


if 6 < 62. The statement of the lemma follows by combining (2.32) and 


(2.36). 
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3. STABILITY RESULTS AT Y = 0. 

Our first statement concerns sufficient conditions for stability 
at zero contact angle: 

Theorem 3.1 ; Let y ~ 0 in (1.1) and assume that the unperturbed 
surface satisfies the Enler-Lagrange equations with the associated 
boundary conditions (2.9), and does not intersect the z-axis. Assume 
further that the function w(s) of (2.3) is twice continuously 
differentiable and that 

w(s) < 0 in some open neighborhoods of s - and s = - (3.1) 

Let K he defined as in (2.15) - (2.17) but with the boundary 
o 

conditions in (2.15) replaced by the fixed end-point conditions (2.19), 

and assume that k > 0. Let d(r\) and d(Z) be defined as in (2.20) and 
o 

(2.21). Then there exists a constant d > 0, such that in (2.4) 

o 

^E(q) > 0 with equality iff T\ = 0 in H (3.2) 

holds for all r\ satisfying the volume constraint = 0 in (2.5), 

the conditions a) - d) of Section 2.2 and the condition 

d(r\) + d(Z) < d . (3.3) 

— o 

Remark : When Y ~ both the energy (1.1) and the liquid volume 

remain unchanged if the liquid-vapor Interface is continued past the 
contact lines by "wetting" dry parts of the container walls. Thus, 
when Y = 0, any configuration is neutrally unstable with respect to 
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such "wetting" perturbations. With our notation, however, the region 
2 is unchanged at "wetting"" (see the beginning of Section 2.2), 
and there is no ambiguity in (3.2) in this respect. 

Proof of Theorem 3.1; Let Q be defined as in (2.10). Then, 
o 

by (2.15) - (2.17) and the representation theorem for quadratic forms 
in Hilbert space (see e.g. Kato [7], p. 322-323): 

inf Q (y) = K >0 (3.4) 

o — o 


where the infimtira is taken over all y satisfying in the conditions 
a) - c) of Section 2.2, (2.11), (2.23) and (2.27). 


Using the notation of (2.7) we put 


f^(n^,s) = f^(ri,s) - f^(^,s) + f (n,s) - Xf^(n,s) 


(3.5) 


where X is the constant in (2.8). Then, by (2.4) and (2.5), 


277 s ref)) 


5E(u) = 


// 


fg(n,s) dsdtf) 


0 s^((f)) 


277 

-// 


fg(w,s) dsdef) 


(3.6) 


for ail n satisfying the volume constraint 6 V(ti) =0. For convenience 
in the following we denote 
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= complement of (EoE^) w.r.t.. E, 


E- = complement of (EnE^) w.r.t. E^ and 


E = EUE = E uE , = EuE_ 
o o + 


and define a function ri(s,(j>) on E by 


n(s,(j)) = 


ri(s,(j)) , (s,<|)) e E 


w(s) , (s,<j)) e , , 


(i.e. Ti is obtained by extending T] by wetting those parts of the 
which dried because of the perturbation r))* By (2.6) and (3.7), 
has the boundaries 


(^) = min|s^((j)) ,sj 


; 0 < ({) < 2tt. 


s^(4>) = max|s^(({)) ,s^[ 


Putting further 


, (s,(j)) e E 


n (s,4>) = 


w 


(s) , (s,({)) e E_j_ , 


v(s,4)) = n(s,(|)) - ri (s,(j>) , (s,$) € E , 


(3.7) 


(3.8) 

wall 

E 


(3.9) 


(3.10) 


(3.11) 


we obtain by (3.6) and (2.5), noting that fg(0^,s) = 0 and f^(0^,s) - 0 
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2ir ((j)) 


6E(n) = 


-If 


0 '^o(4>) 
2tt St (4>) 

=// t 

*'0 •'s (A) 


fg(^,s) dsd<{) 


27T 

■II 


0 •'A(J) 


fg(w,s). dsd<}» 


fp.(n* + V, s) - dsd([) , 


(3.12) 


for all n satisfying the volume constraint 


2ir s,((J)) 


6V(n) 


... 

=J j {fyClI* + V,s) - f^(n*,s)[ dsd(j) = 0. 


(3.13) 


0 


In the sequel and d^, i=l,2,3,..., will denote positive constants, 

independent of u, p, E, s and By (3.8) - (3.11), rj, q*, v and 
il(<})) will satisfy the requirements a) - d) of Section 2.2. Furthermore, 
since w(s) = 0(s - s^) in the neighborhood of s == s^, we may find some 


and d^, such that 


d(Z) + d(n*) + d(n) + d(v) < M^|d(E) + d(n)f 


(3.14) 


if only d(Z) + d(n) < d^. 

By (2.9) and since R(s) >0ins <s<s, the functions R and Z 

o — — 1 

can be continued analytically to some open interval containing [s^,Sj^l 

and it will hold R(s) > R . >0 there. It then follows from (2.7), 

— mxn 

(3.5) and (3.9) that there exist some d 2 and d^, such that f^(ji,s) 
and f^(jl,s) are analytic functions of the arguments n, and in 
the region |r|[ + |ng| + [ri^| ^ ^2 points of E, if only d(E) < d^. 

Hence, putting in (3.12) and (3.13) 
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fg<H* + v,s) - fg(n*,s) 

and 

A 3f„ 

-’®^ “ W (0,s)v + h^(n*,v',s) (3.16) 

there will exist some M 2 and d^, such that 

]hg(Ti*,v,s)i < M2|d(n) + d(E)[ (Vg + v| + v^) 

lh^(H*,V,s)j < M 2 |d(h) + d(E)[ (IVgl + Iv^I + |v]) (3.17) 

holds at all points of S where- jl* and V are continuous, if only 

d(n) + d(Z) < d^. 

We now observe, that since w(s) < 0 in some open neighborhoods of 
s = s^ and s = s the first term to the right in (3.15) gives a ‘ ' 
non-negative contribution to 6E(r|) . To see this, first note that 
= 0 in Z^, by (3.10), (3.5), (2.7) and (2.8). Second, 
when (s,<|)) S we have by (3.10), (2.12) and (2.14) 

fEi]^n*»s)'^v = jaR(s) + Dj^(s)}w’(s)Vg + |c(s) + D(,(s)fw(s)v (3.18) 
where, for some M^ and d^, 

|Dj^(s)[ |Dj^(s)j <M^|jw(s)| +w'(s)^} 



66 


if only d(2) < d^. Furthermore, since w"(s) is continuous in 
will be continuously differentiable there. Hence, after partial 
integration of the first term in (3.18) noting that V(s,<J>)w’(s) = 0 
on the boundaries of we obtain 


27T 5^ ((f)) 




^ ** 

If 

0 

/V 

= / / |~aR(s)w"(s) +D(s)[v dsd(f) 

•'0 «^( j ) 


(3.19) 


where, for some dg, M^, [d(s)| < M^jw’(s)| if only d(E) < dg. Now by 

(3.1), since w(s^) = w’ (s^) = 0, i=0,l, w"(s) < 0 in some open 

neighborhoods of s = s and s = s, , whence the first factor of the 

o 1 

last integrand will be non-negative, if only d(E) is sufficiently 
small. Since further, by (3.8) - (3.11), v(s,(f)) > 0 at interior 
points of Ef> we can then find some d^, such that for all permissible 
r) for which 5V(n) = 0 


2tt 5^((j)) 

<SE(n) > I / f '{A(s)Vg + B(s)v| + C(s)V^ + 2hg(n*,V,s)| 

s ((f))' 
o 


dsd(f> 


(3.20) 


if only d(ri) + d(E) < d^. Furthermore, (3.20) holds with equality 
if E^ is empty or if w(s) = 0 in E+. 

In (3.13) and (3.20) we introduce a change of variable by 


putting 
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s = s(s’,(J)) = s ((j)) +- ^(s- (<J)) - s (<j))) (3.21) 

O S- — S X o 

1 o 

u(s’,(fi) = V(s(s’ . (3.22) 

(3.21) takes Z onto the rectangle in the (s' ,0) -plane. It follows 
by (3.17), (3.21), (3.22) and the smoothness properties of A, B and 
C that in the notation of (2.25) 

<SE(ti) > -|q^ (a + 6A,B + 6B,C + 6C,JJ) 

where 6A(s',4>), 6b(s',<{>) and 6C(s',({») are bounded and integrable 

on Z and such that for some M_ and d^ 
o 5 o 

■|I6a!I^ + IISbII^ + Il5cll^ < M 5 |d(n) + d(Z)[ 

if d(n) + d(Z) < dg. Similarly, by (3.13), (3.16), (3.17) and 
(2.7), there exist Mg and dg, such that 

|(R,y)^I < Mg|d(n) + d(Z)[|[IuJI^ + +llvllj 

if d(n) + d(S) < dg. Furthermore, u(s’,<i>) satisfies (2.23) and 
requirements a) - c) of Section 2.2 on Z^. We may then use (3.4) 
and Lemma 2.2 to conclude that there exists some d^ such that, e.g., 

6E(n) > 3 ^- ilpll^ 

for all n satisfying the volume constraint, if only d(ri) + d(Z) < d^. 

Since, by (3.7) - (3.11) and (3.21) - (3.22). liyil = 0 if and only if 

o 
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n = 0 in Z, this completes the proof of Theorem 3.1. 

The second statement of this section is concerned with necessary 
conditions for stability at zero contact angle, based on the fixed 
end-point conditions (2.19). As may be expected, these will apply 
regardless of the additional ;condltion (3.1); 


Theorem 3.2 ; Let y = 0 and assume that the unperturbed surface 

satisfies (2.9) and does not intersect the z-axis. Let K be defined 

o 

as in (2.15) - (2.17) but with the boundary conditions (2.19) and 

assume that K < 0. Then, for any d > 0 we may find a function n 
o o 

satisfying a) - d) of Section 2.2, the volume constraint &V(T)) - 0 

in (2.5) and the condition d(r\) + <?(ZJ < d , such that in (2.4) 

— o 

5E(n) < 0 . (3.23) 

Proof; We note that the infimum K in (3.4) under the conditions 
o 

stated there is attained for y = ft where y is either y^^^^ or ^1^10 ^2^20 

in the notation of (2.15) - (2.17). We can use y to construct a function 

n with the properties required in the theorem as follows; Let e.g. 

= 2s /3 + s^/3, = s /3 + 2s^/3 and put 

00 o 1 10 o 1 


(s - 


%0 ^ ^^10 


- s) 


00 - 




g(s) = 


(3.24) 


otherwise . 


Then put 


n(s,<|)) = ey(s,4») + a(e)g(s) 


(3.25) 
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where a{£) is chosen so as to make n(s,<{>) satisfy the volume constraint. 

Using (2.5), the assumption (R,1 j)^ = 0 and the contraction mapping 

theorem it can he shown that a(e) is well defined when je| is small 

and that a(e) = o(e), e 0. Since y and g are zero on s = s and 

o 

s = Sj^, as defined in (3.7) will be empty. Hence, by putting 


n*(s,<t)) = 


0 ; (s,4.) e S 

w(s) - ey(s,<f)) ; (s,^) e E 


(3.26) 


we obtain by (3.6) 


2tt s 



1 


U* + a(e)^,s) dsdc]) 


(3.27) 


o 

if only e is small enough for g(s) to be zero within E_, cf, (3.24), 
(3.25). Noting that jr|*(s,(})) | < |eO(s,(f))| in E_ and that the area 
of E_ tends to zero as e 0, we get from (3.27) 


2tt Si 

6e(h) = ey J |u%^^(0,s)y dsd<t> + o(e^) 

0 s ' 


2 ^*^0 


e 0. 


Since K <0, the statement of Theorem 3.2 follows, 
o 

Finally, the following theorem states that if the analytic 
continuation of the equilibrium liquid-vapor interface penetrates 
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the container walls at either of the three-phase contact lines, the 
configuration is unstable : 

Theorem 3.3: Let y = 0 and assume that the unperturbed surface 

satisfies C'2,-9-)- and does not -intersect the z-axis. Assume- further j that 
for i=0 or i=l the function w(s) in (2.2) changes sign at s = s. and 
that w' Hs) is continuous and monotone in some open neighbourhood of 
that point. Then, for any d^ > Oj we may find a function n satisfying 
a) — d) of Section 2.2, the volume constraint bV(r\) - 0 in (2.5) and 
the condition d(vi) + d(Z) < dg j such that in (2.4) 

5E('!\) < 0 . 

Proof: Assume e.g. i=1 . Then, by the mean-value theorem, we may 

find some Sg > s^ such that e.g. 

IWCs)! < (s-s^)v' '(s) 
lw(s) I £ (s-s^ )^w’ ’ (s) 

holds for s^ £ s < Sg . Hence, noting the upper bound for |D(s)| in the 
integrand in (3.19) » for any given d^ > 0 there exists a k > 0 such that 

k 5 3 dQ , (3.28) 

max { lw(s) I , Iw' (s) I } < ^ d^^ < s 5 + k , (3.29) 

-a R('s) w*'(s) + D(s) <0 ; s^<s <s^+k . (3.30) 

Now consider the set of axially symmetr-ie perturbations n = r)(s) for which, 
of. (2.2), 

SQ('j>) = Sq ; s^('Ji) = S^'+ k ; 0 ^ < 2tx . 

Let v(s) be defined as in (3.7)“(3.11) and, keeping k fixed, let d(v) 0 
(i.e. let the perturbation tend to the unperturbed surface extended by 
wetting the portion s^ < s s s^ + k of the container- wall ) . By (3.30) the 
integral (3.19) vill be negative and linear in v while by (3.12), (3.15) 
and (3.-17). other contributions to 6E(n) are, of higher order in v. Hence, 
proceeding as in the proof of Theorem 3.2 to satisfy the volume constraint, 
for any given d.^ > 0 we may choose a function n(s) suitably to satisfy 
6V(n) = 0 and 
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d(v) < , (3.31) 

6E(n) < 0 . 

Since, "by (3.28) and (3.29), using the definitions (2.20), (2.21) and 
(3.7)-(3.11), 

d(Z) = k < I dQ 

lim d(n) = sup 1 ^( 3 ) I + sup iw'(s)| < | d^ , 

d(v)-0 s^i3i3^+k s^^s£s^+k 

the statement of the theorem follows by choosing d^ in (3.31) sufficiently 
small. 

Remark : We note that the instability stated in, Theorem 3.3 is present 

regardless of the stability conditions based on the second variation of (1.1) . 
The method of proof suggests that the instability should show by a thin layer 
of liquid building up past the contact line s = A verification of this 
kind of instability by experiments would of course be virtually impossible 
because of the several idealizing assumptions inherent in the model (cf . 
however [3] where a series of low-gravity experiments performed at the NASA 
Lewis Research Center are studied computationally) . 
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ABSTRACT 

•We study computationally the stability, under gravitational and 
surface forces, of a liquid in a circular cylindrical container with a 
concave spheroidal bottom, for the case in which the volume of liquid 
is sufficiently small so that the bottom is not covered entirely. We 
assume the gravitational field to be directed along the axis of symmetry 
of the container, and for a specific container shape we compute the 
critical Bond number as a function of liquid volume for contact angles 
Y = 0°, 1°, 2°, and 4°. For the case y = 0® we present graphically 
several critical equilibrium configurations and corresponding pertur- 
bation modes. 
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1. STATEMENT OF THE PROBLEM 

In this paper we present the results of a computational study of 
the stability of a liquid in a rotationally symmetric container subject 
to gravitational and surface forces. We consider vertical right cir- 
cular cylindrical containers with concave spheroidal bottoms, for the 
case in which the volume of liquid is sufficiently small so that liquid 
lies only in an annular region of the container C^ig. 1) • We are 
interested specifically in the case for which the contact angle y is 
zero, or nearly zero, and our numerical results are obtained for a con- 
tainer currently used for the storage of liquid fuels in National 
Aeronautics and Space Administration Centaur space vehicles, for which 
the axial ratio of the bottom is b/a = 0.724. 

A vertical section through the axis of the container is depicted 
in Fig. 1, along with the associated cylindrical coordinate system. 

The container may be in motion, but the net external gravitational 
force is assumed to be uniform and directed parallel to the axis of 
symmetry. It is well known, that even if the gravitational force is 
directed upward, liquid may be in stable equilibrium at the container 
bottom because of the effect of surface forces. For a given liquid 
voliime, stable configurations of this kind are possible only if the 
magnitude of the upward-directed gravitational force does not exceed a 
certain critical value. This critical value depends on physical para- 
meters such as the liquid- vapor surface tension coefficient, the differ- 
ence in liquid and vapor densities, the liquid- container contact angle. 
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Figure 1 : The container and the associated cylindrical 

coordinate system. The bottom has the shape 
of an ellipsoid of revolution with axial ratio 
b/a. 
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and geometrical parameters such as the container size and shape. The 
combined effect of certain of these parameters may be represented by 
the dimensionless Bond number (see (2.11) below), whose critical value 
for our problem is a function only of the container shape, the 
liquid volume, and the contact angle. In this study we determine 
con 5 >utationally the critical Bond number as a function of the liquid 
volume for fixed contact angle and container shape. 

Our approach is that of static analysis , i.e. we consider the 
total potential energy of the liquid-container system (in a container- 
fixed frame of reference) , given by 

E = a(A^ - cosy) + Eg (1.1) 

(cf. Reynolds and Satterlee [3], p. 394-396). Here a > 0 (the liquid- 
vapor surface tension coefficient) and 0 < y < it (the contact angle 
between the liquid vapor interface and the container wall and bottom) 
are constants determined by physical properties of the liquid and the 
container, A^ and A^ are the areas of the liquid-vapor and the liquid- 
container interfaces, respectively, and E is the gravitational poten- 
tial energy of the liquid. A configuration of liquid is in stable 
equilibrium if and only if the total potential energy (1.1) is minimal 
compared with that of any nearby configuration having the same liquid 
volume . Thus the critical Bond number for a certain volume of liquid 
is the one at which E in (1.1) ceases to have a strict local minimum 
with respect to all perturbations that conserve the liquid volume. In 
Section 2 below we give a summary of an analysis of this constrained 
minimization problem using methods of variational calculus. 
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Figure 2: The arc-length, normal-displacement coordinate system. 
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continuation o£ the unperturbed surface to the left in Fig. 2 . 

(The continuation may be done in several ways, perhaps most conveniently 
by using the fact that RCs) and Z(s) will be analytic- functions for 
equilibrium configurations, cf. (2.9] below.) 

In a way similar to (2.2), the container wall and bottom are 
described by 


r = R(s) - w(s)Z* (s) 


0 < (j) 2 tt 


(2.3) 


' z = Z(s) + w(s)R'(s) 

in some neighborhoods of the unperturbed contact lines s = Sq and s = s^^. 
Then clearly w(s^) = w(s^) = 0. The function w(s) will depend on the 
shape of the container wall and bottom and, in 5 )licitly, the shape of 
the unperturbed liquid-vapor interface. The representation (2 .'3) is 
convenient for the purpose of deriving the differential equations; 
however, in the actual computations we must, of course, make use of the 
known, configuration- independent shape of the container, cf. (3.4) and 
(3.5) below. 

The increments of the total potential energy E and the liquid 
volume V caused by the perturbation r)(s,4>) in (2.2) may then be computed. 
We obtain in a straightforward way 
<SE (n) = a (SA^ - SA^ cosy) + SE^ 


27T s ((|)) 


0 Sq((1)) 


2ir 


0(COSY • fj 5 ^(w,s) - f^(0,s)) + pg £g(w,s)| ds d(j> ^ 

0 A . 


<i> 


(2.4) 
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2ir 2ir 

6VCnD = /* f fwCa^s) ds d(|) - f f f„Cw,s) ds d(|) , 

0 Sq(c|)) 0 (2.5) 

where we have put 

T 

n = nCs,4>D = CnCs,4>), 'n^Cs>4>D, n^(s,cf))) , 

w = w(s) = (w(s), w'Cs), 0)^ C2.6) 

= the interval CSqC4>D> s^) U s^C4>D) , 

and where the functions fyCas^Dj and £^( 11 , 5 ) are given by 

(denoting R = R(s) , Z = Z(s), and C = Z' 'R' -R''Z' for brevity) 

f^Ca.s) = { (R - nZ’)^Cri3 + (1 - nc)^) + n^Ci - nC)^ } " 
fyCn>s) = pjR - |(z' + RC)n + yZ'Cn^} 

(2.7) 

fg(n,s) = n{RZ + ^(RR' - zz> - Rzc)n - 

- -|((RR' - ZZ')C + R'Z’)n^ + jR’Z'Cn^} 

In (2.4) p, the liquid density (or, more precisely, the difference 
between the liquid and vapor densities) , is assumed to be constant, and 
g is the gravitational constant (which may assume any value) defined so 
that g > 0 if the gravitational force pulls toward the negative z-axis. 
The condition, that all first-order p-terms in 6E(p) as given by 

(2.4) should vanish for all p such that 6V(p) = 0 in (2.5) is then 

3f, 9£ 3£„ 

^C0,s) + pg 3^(0,s) - X 3^(0,s) = 0 

in Sq < s < Sj , with the side-conditions 

cosy • £^(w(sp,s^) - £^(0,s^) = 0, i = 0, 1, 


( 2 . 8 ) 
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where A is a constant C^he Lagrange multiplier) . Putting here 
B = Pg/cf, H = X/ a, and using C2-7) and the identity 

R» (s)^ + Z' Cs)2 = 1 


(which holds because s is the arc-length) (2.8) becomes the Euler - 
Lagrange boundary value problem for the equilibrium liquid-vapor 
interface 

R' ' = -Z' (BZ - H - Z'/R) 

Z'» = R’ (BZ - H - Z'/R) 

(2.9) 

w(Sq) = w(Sj^) = 0, 

w’ (Sp) = -tany, w* (s^^) = tany 


In general, the requirement that (2.9) should have a solution for a 
given container shape will restrict the B, H, and V - values to some 
two-dimensional subset of the (B,H,V) - space. 

The equations (2.9) have the following invariance property : if 

A is any positive constant and if we put 

R(s) = |;R(As), Z(s) = |z(as),w(s)= ~w(j!.s) (2.10) 

then R and Z will satisfy (2.9) with B replaced by B = , H 

replaced by H = H5, and w replaced by w . The transformation (2.10) 
means simply that R , Z , and w describe a liquid -tank configuration 
obtained by uniformly enlarging the original one by a factor l/l. 
Therefore, if we define dimensionless constants B^ and by 


"o' 




a 


«o ' » " = # 


( 2 . 11 ) 
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where a is the container radius, (the Bond number) and will be 
invariant under uniform re-scalings of equilibrium liquid-tank config- 
urations . We will therefore present the results of our coji^jutations 
below in .terms of the B^ in (2.11) in order to facilitate their use for 
arbitrary- sized containers. 

We assume now that the liquid- vapor interface (2.1) satisfies the 
Euler- Lagrange equations (2.9). Then the condition, that all second 
order q-terms give a non -negative contribution to 6E(q) in (2.4) for 
all n such that in (2.5) 6V(q) = 0, may be written as 


2tt 

//I 




0 s. 


2ir 


J { •+ aj^n(s^,({))^ I 


d<J) > Q 


0 

for all q. = q(s,(j)) .such that 

27T s. 


// 


R(s) q(s,({)) ds d<fi = 0 


0 s. 


( 2 . 12 ) 


(2.13) 


In (2.12) we have denoted 

A(s) = -2R*' + b|rR' - ZZ* + ZR(R"Z' - Z*’R')j- 
+ h|z» - R(-R* 'Z' - Z' 'R’)J 

and 

a = ~ ) cosy * f.(w(s),s) - f.C0,s) + Bf (w(s),-s) 

^ tan^Y da ( ^ ® 


(2.14) 


- H£^(w(s),s) 


i = 0, 1 


(2.15) 


s=s. 
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It can be shown that if (2.12) (with the side-condition (2.13)) holds 
with strict inequality for all nonzero functions q(s,4>) then in fact 
E will be locally minimal at q = 0, i.e. the configuration described by 
the R(s) and Z(s) in (2.9) is a stable one. It can further be sho\m 
that if (2.12) (under (2.13)) does not hold for all q, then E' cannot 
be locally minimal at q = 0 and the configuration is an unstable one. 

Thus the critical value of B will be the value at which the transition 
between these two cases occurs (provided that the corresponding solu- 
tions R(s), Z(s) to (2.9) describe an equilibriinn liquid- vapor interface 
that is physically realizable, which may not always be the case - see 
the end of this Section). When B is critical in this sense, (2.12) 
still holds (vmder (2.13)) but there exists some non-zero q = q(s,(|)) for 
which (2.12) holds with equality. 

The inequality (2. 12)- (2. 13) may be analysed in terms of an 

associated sequence of eigenfunctions of the form 

00 

{p.,(s) cosk(j) } (2.16) 

i=l,k=0 

CO 

with eigenvalues given by the solutions to the 

eigenvalue problems 

I - ^ {RCs)y^^(s)} - * ACs)}- y,^Cs) = 

RCs.Jy' Cs.) = C-l)^a M,, (s.) j - 0, 1 

3 ^ 3 1 k = 0, 1, 2, ... 

i = 1, 2, 3, ... 




C2.17) 
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where A(s), and were introduced in (2.14)-(2.15) . It can be shown 
from (2.17} that the eigenvalues increase with the index k. It further 
follows from (2.16} that all eigenfunctions but those with k = 0 in (2.16) 
satisfy trivially the condition (2.13). Then it follows, that (2.12} 
under the side-condition (2.15} holds if and only if (assuming to be 
ordered increasingly with i} 

min|x,,,E BJX^J > 0 , (2.18} 

where, because of (2.13), solutions to: 

00 2 

Minimize E 3-X. under the constraints 
■ 1 - 10 . 

1=1 

oo 2 . ■ 

V g. = 1 and (2.19) 

i=i ^ 

CO 

E >R) = 0 . 

i=l ^ 

In (2.19) we have denoted 

2ir Sj 

(f,g) = // f(s,4>) g(s,«})} ds d4> , 

f 

and assumed (u. }. , to be normalized to (u. ,u. } = 1, i= 1,2,3,... . 

In all the cases studied below it will in fact hold | (ij^q,R) | » 
|(y^p,R}j, i=.2,3,... and CX 2 Q - X^^q) » C^^l ~ ^10^ ^ so- that the minimum 
in (2.18) is attained at X^^^. Thus, the critical value of B is determined 
by the condition X^^^ = 0 (with the above mentioned restriction concerning 
non-realizable equilibriiim configurations). By (2.14), (2.15), and 
(2.17), we obtain that this condition is equivalent to 
requiring that the Jacobi -Legendre boundary value problem 
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- Ry'-'Cs) - R'li'Cs) +1^- 2R" + B{RR' - ZZ' + ZRCR'*Z> - Z''R')) + 

+ H{Z* - RCR'»Z' - Z’'R*)}1 u(s) = 0, s. < s < 

^ ^ ^ (2.203 

(-1)^ tany-y' (s^) = {sin^Y* (R' ’2' ~ Z'’R’) - cos^Y'W*} *y(s.) 

s=s. 

1 

i = 0, 1, 

( 2 . 21 ) 

should have a non-trivial solution y(s) . 

We remark that the formulas (2.12) and (2.15) in the above dis- 
cussion are meaningful only if the contact angle y is strictly positive. 
However, one of our principal interests is .the limiting case y = 0. 
Therefore a special analysis is needed in order to determine the proper 
limiting form of the above conditions (2.19)- (2. 20) when y ^ 0. It can 
be shown (Karasalo* [4]) that if w'*(Sj.) < 0, i = 0, 1, then the dif- 
ferential equation (2.20) with the fixed end-point boundary conditions 

U(Sq) = y(sp = 0 (2.22) 

is the correct one to use when y = 0. Furthermore it holds when y = 0 
that ^ in the set of solutions R(s), Z(s) to (2.9) obtained by keeping 
the volume fixed and varying B (and H, cf. the comment after (2.9) above) 
w* * (Sq) or w* ’ (Sj^) change sign as functions of B at some value of B , 
then this B- value is critical even if (2. 20) -(2. 22) lacks nontrivial 
solutions . This is so because only solutions R(s), Z(s) for which 
•w' * (Sp) < 0 and w’ ' (s^^) < 0 hold are permissible for y = 0 due to the 


constraints imposed by the container geometry. Our computations show, 
in fact, for the Centaur space vehicle example, that small liquid 
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volumes become unstable because of the conditions (2,20)- (2.22 j , 
whereas the stability of large liquid volumes is decided by the con- 
straint w* * (Sj^) < 0 . The transition between these two conditions occurs 
at a certain well defined volume, c£. (3,6)- C3. 7) below. 

5. COMPUTATIONAL PROCEDURE 

With a given liquid volume in a given container of the shape 
shown in Fig. 1 we associate a dimensionless fill height , defined as 
follows: let be such that the given volume V coincides with the 

volume bounded by the container wall and bottom and the plane z = . 

Then the fill height for the volume V in the container with 

radius a is 



We shall compute the critical Bond number (cf.(2.11)) as a 
function of h^. Obviously, by (2.10)- (2.11) and (3.1), these quantities 
are invariant under uniform re-scalings of the container, and we can 
therefore restrict our computations to a container with a specific 
radius, e.g. a = 1. We are interested only in the fill-height range 
0 < hy < b/a, i.e. only in volumes that are smaller than that of the 
annular crevice at the container bottom. 

- Before describing our computational algorithm in detail we shall 
give the explicit form of the boundary conditions (2.21) at s = s^ and 
s = Sj, respectively. In a neighborhood of s = s^ there holds by 
(2.3), since the bottom is an ellispoid of revolution (cf. Figs.l and 2) 
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CRCs) - w(s)Z»Cs))^ CZCs) + wCs)R«rs))^ , 

2 ,2 - 1 . 

a b 

Similarly, in a neighborhood of s = we have 

RCs) - w(s)Z‘(s) = a 

By differentiating these expressions twice and using (2.9), (2.21) 
becomes after some straightforward manipulations (we put 
= R(Sj^)j i = 0, 1, for brevity) 


(3.2) 


(3.3) 


sinyy' (s^) 


f a"^b"^ 


cosy (BZq 



•y(So) 

(3.4) 


siny-y»(s^) = -cosy (BZ^^ - H - ) d(Sj^) . (3.5) 

Our computations axe carried out for the case b/a = 0.724 and 
proceed in two principal steps. In the first of these (which requires 
the main part of the computational effort) we determine successively 
some 50-60 points on the curve non-equidistant values 

of hy. Each of these points is obtained in the following manner: 

a) We choose a fixed point R = a = 1, Z = Z^ at the cylindrical 

container wall and ’’guess'* a corresponding pair of values for and 
Hq in a way to be specified below. . 

b] ' We put R(Sj^) = a = 1, ^(sp = 2^, R’ (s^^) = siny, Z’(Sj^) = cosy 
to satisfy the boundary conditions at s = in (2.9) (we may choose 

arbitrarily, e.g. s^^ = 1) . We further choose a pair of values 
y’(Sj), v(Sj), not both zero, consistent with the boundary condition 
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(3.5) (except for the large cases for Y = cf. (3.6) below). 

c) We solve simultaneously the differential equations (2,9) and (2.20) 
numerically, integrating from s = s^^ backwards with a standard fourth- 
order Runge-Kutta scheme. The stepsize of the integration is kept 
constant except for the last step which is adjusted (using a secant 
method) so as to make the last computed point of the solution to (2.9) 
lie on the container bottom profile (for the case y = 0 we adjust the 
last step so as to make the normal of the computed solution to (2.9) 

at s = Sq intersect the bottom profile at an angle of 17/2) . Thus we 
have ensured w(Sq) = 0 (w' (s^) = O.in the case y = 0) . 

d) We compute the discrepancies in the boundary condition (3.4) and 
the remaining boundary condition at s = s^ in (2.9). We adjust Bq and 
Hq (using eventually a Newton-type method to obtain the corrections) and 
repeat from b) above until the corrections in and Hq are less than 

a prescribed tolerance. 

e) We make a final integration computing this time also the liquid 
volume V, simply by adding the appropriate extra differential equation 
to the others. Then the corresponding h^- value is obtained using 
(3.1). 

We repeat -the steps a)-e), using a set of some 50-60 regularly spaced 
Z^-values. To obtain the initial "guesses" for Bq and Hq in a), we 
extrapolate the functions Bq^ = Bq^(Zj^) and Hq^ = Hq^(Z^) to the next 
Zj^-value, fitting two quadratic polynomials in Z^ through the three 
closest previously computed values of log|BQ^| and log[HQ^|, respectively. 
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The functions log|BQ^(Z^)j and loglHQ^CZ^^^ | turn out to be close to 
linear in Z^, and the accuracy in the guessed values was found to be 
very good fthe "guesses" have in general 3-4 correct decimals when the 
spacing of the Z^^-values is 0.02a). 

For all the contact angles studied, the computed points on 
the curve = ®0c^^V^ indicate that logjBp^l is only mildly nonlinear 

as a function of h^. The second main step of our computation is to 

fit a cubic spline through the computed points on the curve logiBg^l = 
log I Bq^ ( hy) I , In this way we obtain a convenient and satisfactorily 
accurate representation of the sought function B^^ = ®0c^^V^ throughout 
the entire h^- interval of interest. 

We have studied the contact angle values y = 4®, 2°, 1°, and O'*, 
the last of these values being the case of main interest. In the case 
Y = 0 we find that the above algorithm must be modified in the following 
way : IWien, in step a) above, Zj^ > Z* (corresponding to h^ > h^, where 

Z* = 0.7014 and h* = 0.5031 with four correct decimals), then the 

condition 

w'»(s^)<0 (3.6) 

(which expresses the condition that the equilibrium liquid-vapor 
interface must lie inside the container for s < s^^, cf. the discussion 
at the end of Section 2) places a more restrictive bound on B^ than 
the conditions under d) above. By (3.3), (3.6) is equivalent to 

R"(sp = -(BqZj, - Hq - 1) < 0 


(3.7) 
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Hence, when y = 0 and > Z| we solve, in steps b) and c) above, only 
.the .differential equations (2,..9.), and adjust, in step d) , and so as 
to satisfy w(s^) = o and - Hg - 1 = 0 (using again a Newton-like 

method to obtain the corrections) . 

4 ■ NUMERICAL RESULTS 

Table 1 shows the cubic spline approximations to the functions 

Bqc = contact angle values y = 2°, and 4° (for 

the case b/a = 0.724). The relative error in each entry shown is less 
-4 

than 10 , as estimated from repeated computations with different step- 

sizes in the numerical integration and different spacings for the 
Z^-values used in step a) in the computational algorithm (cf. Section 3) 

Figures 3-6 show graphs of the functions in Table 1. For practical 
reasons we use a logarithmic scale on the Bg^-axis, for which the curves 
are close to linear. Figure 7, showing all the graphs from Figures 3-6 
simultaneously, illustrates the almost insignificant dependence on y in 
this range. 

In Figure 8 we show the equilibrium surfaces at critical Bond 
number for the fill heights hy = 0.2(0. 1)0. 7 for the case y = 0°. The 
curves were plotted by numerical integration from starting points at 
the cylindrical wall using the subroutine package GRAFPAC available at 
LBL for graphical display of the results. The starting points for the 
integrations were obtained using cubic spline fitting to the points on 
the curve Z^ = Z^^ (hy) , which are known as a "by product" of the 
Boc = ®0c^^V^ ~ calculation (cf. steps a) - e) of the algorithm described 
in Section 3) . 
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The graphs of Figure 9 represent solutions to the Jacohi- Legendre 
equation (2.20) in the cases corresponding to the equilibrium configura- 
tions shown in Figure 8. The functions shown in Figure 9 are of the 
form 

where y(s) is a solution to (2,20) satisfying p(Sq) = 0. The function 
v(s) depicts the radial dependence of the perturbation p(s) cos(J>, but 
in terms of vertical displacement (whereas y(s) is the radial pertur- 
bation profile in terms of displacement normal to the unperturbed 
surface , cf. (2.2)). The abscissa of the graphs in Figure 9 is R(s) 
and the functions are normalized so that dv/dR » 1 at the left end 
points. Theoretically, for h^ >h* 0.5031 (see (3.6)) the functions 

V(s) defined in (4.1) have a singularity at s = s^, i.e. , at R = 1, 
while for hy < h* v(Sj^) is nonzero but finite. Hence the graph cor- 
responding to hy = 0.5 in Figure 9 in fact intersects the line R = 1, 
whereas the two graphs above it do not. 

Finally, in Figures 10-15 we show the equilibrium liquid-vapor 
interfaces of Figure 8. (solid lines) together with the equilibrium 
liquid-vapor interfaces superimposed by a small multiple of the cor- 
responding v(s) as given by (4.1), These curves are of some interest 
because, theoretically, when Bq = and hy < h* the only possible 
initial shape of an unstable perturbation is given by v(s) cos(j> (in 
terms of vertical displacement) . For hy > h^ it can in fact be shown 
(Karasaio [5]} that the configurations are unstable at for. 
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e.g.j perturbations that build up some suitably chosen, axially sym- 
metric layer o£ liquid above the unperturbed contact line at the cylin- 
drical wall. Notice, however, that the dashed curves in Figures 14 and 
15 do not show such a perturbation, but they are nevertheless included 
here for completeness. These results on the initial perturbation shape 
rely, of course, on several idealizing assumptions, such as that (1.1) 
holds exactly, that it represents all boundary constraints, that 
viscosity effects need not be included, etc., and the conclusions from 
Figures 10-15 should not be drawn too far. 

The numerical results presented here have been found to be consist- 
ent with preliminary experimental results obtained at the NASA Lewis 
Research Center [6] . 
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Table 1. 

Critical Bond 

mimber, , as 

Oc 

function of fill height, h„, at 


Y = 1% 2 

% and 4“. 




II 

O. 

o 

Y = 1 ° 

Y = 2“' 

- Y = 4“ 

0:200 

-480.4283 

-480.3526 

-480.1183 

-479.1676 

0.210 

-412.6934 

-412.6251 

-412.4212 

-411.5999 

0.220 

-356.8787 

-356.8212 

-356.6472 

-355.9423 

0.230 

-310.4444 

-310.3955 

-310.2463 

-309.6393 

0.240 

-271.5210 

-271,4787 

-271.3503 

-270.8253 

0.250 

-238 ..6539 

-238.6178 

-238.5075 

-238.0536 

0.260 

-210.7102 

-210.6795 

-210.5850 

-210.1928 

0.270 

-186.8042 

-186.7783 

-186.6978 

-186.3599 

0.280 

-166.2339 

-166.2123 

-166.1443 

-165.8546 

0.290 

-148.4377 

-148.4200 

-148.3635 

-148,1170 

0.300 

-132.9638 

-132.9499 

-132.9038 

-132.6967 

0.310 

-119.4456 

-119.4352 

-119.3990 

-119.2280 

0.320 

-107.5834 

-107.5765 

-107.5497 

-107.4123 

0.330 

-97.1310 

-97.1274 

-97.1098 

-97.0041 

0.340 

-87.8844 

-87.8843 

-87.8757 

-87.8004 

0.350 

-79.6741 

-79.6776 ' 

-79.6781 

-79.6321 

0.360 

-72.3581 

-72.3654 

-72.3752 

-72.3581 

0.370 

-65.8171 

-65.8287 

-65.8482 

-65.8598 

0.380 

-59.9504 

-59.9666 

-59.9965 

-60.0369 

0.390 

-54.6722 

-54.6937 

-54.7349 

-54.8046 

0.400 

-49.9096 

-49.9374 

-49.9908 

-50.0906 

0.410 

-45.6000 

-45,6353 

-45.7024 

-45.8332 

0.420 

-41.6896 

-41.7341 

-41.8166 

-41.9796 

0.430 

-38.1320 

-38.1879 

-38.2878 

-38.4843 

0.440 

-34.8869 

-34.9574 

-35.0769 

-35.3082 

0,450 

-31.9190 

-32.0084 

-32.1502 

-32.4173 

0.460 

-29.1976 

-29.3120 

-29.4786 

-29.7822 

0.470 

-26.6956 

-26.8433 

-27.0371 

-27.3775 

0.480 

-24.3887 

-24.5813 

-24.8043 

-25.1807 

0.490 

-22.2554 

-22.5086 

-22.7615 

-23.1725 

0.500 

-20.2759 

-20.6107 

-20.8925 

-21.3355 

0.510 

-18.4509 

-18.8751 

-19.1830 

-19.6544 

0.520 

-16.8150 

-17.2908 

-17.6204 

-18.1156 

0.530 

-15.3484 

-15.8474 

-16.1930 

-16.7068 

0.540 

-14.0314 

-14.5345 

-14.8900 

-15.4170 

0.550 

-12.8467 

-13.3417 

-13.7013 

-14.2359 

0.560 

-11.7792 

-12.2587 

-12.6174 

-13.1545 

0.570 

-10.8157 

-11.2756 

-11.6291 

-12.1641 

0.580 

-9.9447 

-10.3828 

-10.7281 

-11.2570 

0.590 

-9.1561 

-9.5717 

-9.9065 

-10.4259 
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Table 1. (cont.) 


0,600 

- 8.4411 

0.610 

- 7.7920 

0.620 

- 7.2017 

0.630 

- 6.6644 

0.640 

- 6.1746 

0.650 

- 5.7276 

0.660 

- 5,3192 

0.670 

- 4.9456 

0.680 ■ 

- 4,6035 

0.690 

- 4,2898 

0.700 

- 4 ,. 0020 


- 8.8342 

- 9.1571 

- 8.1630 

- 8.4730 

- 7. 5517 

- 7.8483 

- 6.9942 

- 7,2773 

- 6. 4854 

- 6.7551 

- 6.0204 

- 6.2770 

- 5.5951 

- 5.8390 

- 5.2057 

- 5.4373 

- 4.8487 

- 5.0686 

- 4.5212 

- 4.7299 

- 4.2204 

- - 4.4185 


- 9.6644 

- 8.9662 

- 8.3259 

- 7.7384 

- 7.1990 

- 6.7035 

- 6.2481 

- 5.8291 

- 5.4435 

- 5.0884 

- 4.7610 
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Figure- 8 ; Equilibrium surfaces at critical Bond number for y 
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Figure 9 : Radial dependence of first angular perturbation 

inodes in terms of vertical displacement. Curves 
normalized so that ~ 0 and ('^)„ = 1. 

O O.K. 

(R =' lef t end point of curve.) 
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Figure 10 : Solid Curve : Equilibrium surface at critical Bond number 

for Y = O'" and = 0.2. 

Dashed Curve: Solid curve superimposed with the corre- 


sponding perturbant mode from Figure 9. 
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XBL 773-7881 

Figure 11 : Solid Curve : Equilibrium surface at critical Bond number 

for Y = 0® and = 0.3. 

Dashed Curve ; Solid curve superimposed with .the corre- 
sponding perturbation mode from Figure 9. 





spending perturbation inode from Figure 9. 
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Solid' Curve : Equilibrium surface at critical Bond number 

for Y = 0® and = 0.5, 

Dashed Curve : Solid curve superimposed with the corre- 

sponding perturbation mode from Figure 9. 
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Figure 14 : Solid Curve : Equilibrium surface at critical Bond number 

for Y = 0° and = 0.6. 

Dashed Curve : Solid curve superimposed with the corre- 

sponding perturbation mode from Figure 9. 
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XBL 773-7878 

Figure 15 : Solid Curve : Equilibrium surface at critical Bond number 

for Y = 0° and = 0.7, 

Dashed Curve : Solid curve superin^osed with the corre- 

sponding perturbation mode from Figure 9. 
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wsm NOT PMliS 


ABSTRACT 


In this paper we calculate the small-amplitude periodic sloshing 

modes of a liquid in a vertical right circular cylinder with a concave 

spheroidal bottom, for the case in which there is not sufficient liquid 

to cover the bottom entirely. Equilibrium free surfaces of the liquid 

were calculated by the program CAPIL for the case in which the ratio 

of the minor and major semi-axes of the spheroidal bottom was 0.724. 

Perturbations about these surfaces were calculated by the program 

SLOSH. For the fill heights that were studied, and to the accuracy 

of these calculations, we found the same critical Bond number, B , 

crxt 

for instability of the free surface as was found in the static analysis 
of P. Concus and I. Karasalo for the same test problem. Furthermore, in 
agreement with their calculation we also foxind no equilibrium surfaces 
for this problem for fill heights greater than 0.503 and for Bond nim- 
bers B < 0 ♦ B"or fill heights ranging from 0.20 to 0.45 we 

found unstable equilibrium surfaces for a range of Bond numbers, 

^conv ^ ® ^crit * frequencies or growth rates were calculated for 
numerous equilibrium surfaces. Growth rates of the maximally unstable 
modes were calculated for fill height 0.30 and various Bond numbers. 
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1. Introduction 

In this paper we calculate the small-amplitude, periodic sloshing 
modes of a liquid in a rotationally symmetric cylindrical container 
under the effect of surface and gravitational forces. We consider 
a right circular cylinder, oriented vertically, with a concave spher- 
oidal bottom, for the case in which there is not sufficient liquid to 
cover the bottom entirely. This is the same configuration for which 
a stability study was carried out in [1] . Numerical results are 
obtained for a container currently used for the storage of liquid fuels 
in National Aeronautics and Space Administration Centaur space vehicles, 
for which the axial ratio of the bottom is b/a=0.724 . A vertical 
cross section of the cylinder and liquid is shown in Figure 1. 

Equations describing the sloshing motion of liquids in rotation- 
ally symmetric containers are derived in [2] using a surface-normal 
polar coordinate system particularly suited to such problems. It is 
assumed there that the fluid flow is irrotational and incompressible and 
the free-surface boundary conditions are obtained from the time- 
dependent Bernoulli equation and the kinematic equation. The dif- 
ference in pressure across the free surface at any point, due to the 
interfacial surface tension, is proportional to the mean curvature at 
that point. The edges of the surface satisfy time- independent contact 
angle conditions with the container bottom and the cylinder wall. We 
follow the derivation in 12] for obtaining the equations of motion for 
the case studied here, but we use a different technique for obtaining 
the numerical solution. 
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Figure 1. Vertical cross section of the cylinder and liquid. 


XBL 781-103 


116 


2. Scaled Variables 

We consider a circular cylindrical coordinate system with the z 
axis along the cylinder’s axis of symmetry. It is convenient to 
define scaled length and" time variaBles . Let symhols with a bar" over' 
them denote the corresponding physical, unsealed variables. Let 
r = r/a 
z = z/a 

t = i [(1 +]b|) CT/pa^]^''^ 

H = H a 

2 2 

B = Ka = pga /a 

where a , the cylinder's radius, is the characteristic length used for 
scaling, t is the time, p is the difference in densities between the 
liquid and gas phases, g is the acceleration due to gravity, consid- 
ered positive when directed vertically downward, Cf is the gas-liquid 
surface tension, K is the capillary constant, B is the Bond number, 

H is the mean curvature at a point on the free surface, considered 
negative when the surface is concave upward, p is the gas pressure, 

g 

and p^~ is the liquid static pressure at the height z = 0 . 

The difference in pressure' across the free surface satisfies the 
equation 

p - Pg = - 2H Cj/a . 

The liquid static pressure is given by 
P - Pq = - Pgza • 

From these equations it follows that the curvature H at any point on 

the equilibrium free surface is related to H , B , and z bv 

o 

2H = 2H + Bz . 
o 
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3. Equilibrium Free Surface 

We consider the vertical cross section through the axis of the 
cylinder shown in Figure 2. The cross section of the liquid is 
bounded by three curves: the meridians along the free surface, the 

cylinder wall, and the container bottom. Let s be the arc length 
along this boundary, increasing clockwise. Let s = 0 be the inter- 
section of the meridians on the free surface and the bottom, and let 
s = S be' the intersection of the meridians on the free surface and the , 
cylinder wall. 

The equilibrium free surface is rotatlonally symmetric about the 
axis of the cylinder. Its height is a function of r only and not of 
0 . Thus the equilibrium surface can be described parametrically by 
the equations 

r = R(s) and z = Z(s) 

for 0 < s < S and 0 < 0 < 2 tt . Let Jp be the angle in the cross- 
sectional plane between the tangent at a point on the free surface and 
the horizontal. Let be positive when the surface slopes upward in 
the direction of increasing s . Then 
tan = Z /R , 

where the subscript s denotes d/ds . Let the spheroidal bottom be 
described by 

" = Zs(r>- 

for R(0) ^ r < 1 and 0 ^ 0 < 2 tt . Let X denote the angle in the 
cross-sectional plane between the tangent at a point on the bottom and 
the horizontal. Let x 1*® negative when the bottom slopes downward in 
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Figure 2. Vertical cross section of the liquid showing coordinates. 
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the direction of increasing r , as it does in our case. Then 
tan X = Zg(r) 

The equilibrium free surface is- the solution of the time indepen- 
dent Bernoulli equation, 

ij; = 2H + BZ - (sin i(<)/R , (3.1) 

s o 

with 


R = cos tI) 
s ^ 


Zg = sin If), 


subject to the contact conditions. 


Z(s) = Z (R(s)) at s = 0 

D. 


R(s) = 1 


at s = S 


. (3.2) 

(3.3) 

(3.4) 

(3.5) 


(the scaled radius of the cylinder is 1), and subject to the contact 
angle conditions, 

ifi”X = Y3^ts = 0 (3.6) 

tt/ 2 - )j) = y at s = S , (3.7) 

where y is the contact angle. The volume of the liquid in the 
cylinder is 


V = 27T 


[Z(s) - Zg(R(s))]R(s) cos if)(s) ds . 

J 

0 


(3.8) 


This last equation determines implicitly the value of of V is given. 

Equations (3.1) - (3.8) are the equations for the equilibrium free 
surface. The solution of these equations varies with the volume, the 
Bond number, the contact angle, and the shape of the bottom of the 
container. Depending on the values of these parameters, ' there may be 
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no, one, or more solutions of these equations [l]; If the equilibrium 
surface exists, it may be stable or unstable to small perturbations. 

These equations are solved by the program CAPIL 1^3] . This program 
uses PASVA2 [4], a general-purpose finite difference solver for non- 
linear first-order systems of differential equations subject to two- 
point boundary conditions. PASVA2 solves these equations by iterating 
from an initial approximation to the surface. Either the user can 
supply the initial approximation, or the subroutine CYLCUR can generate 
it. When making calculations with the same fill for a sequence of Bond 
numbers, we let CYLCUR generate the initial approximation for the first 
case and use the output of each case as the initial approximation for 
the next case. 
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4. Small-Amplitude Periodic Sloshing Modes of the Liquid 

The sloshing motion is treated as potential flow in an incompres- 
sible fluid. The fluid velocity v at any point is the gradient of 
a potential function (}) 

V = V<) . 

Since the fluid is incompressible, V *v = 0 , so <|) satisfies Laplace’s 
equation 

At'= 0 . (4.1) 

The boundary condition on 4> along the cylinder wall and ,the bottom 
is 


= 0 , (4.2) 

where the subscript n denotes the outward normal derivative. 

The displacement of the free surface from its equilibrium will be 
described in surface polar normal coordinates s , 0 and n . The 

coordinate s is the arc length along the equilibrium surface, and 
the coordinate ri is the displacement normal to this surface [l] . The 
perturbed surface is described by 
n = H(s,0,t) . 

The time-dependent Bernoulli equation is linearized in the perturbation 
H . Since is the component of fluid velocity normal to the equi- 

librium surface, H and ^ are related by the kinematic equation on 
the eqiulibrium surface 



This is the boundary condition on <f> 
on the unknown function H . 


(4.3) 

along the free surface; it depends 
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The sloshing motion will be analyzed in terms of normal modes 
^ = ^(r,z) cos(m0) cos(o)t) , 

H = H(s) cos(m6) sin(wt) . 

Equation (4.3) can be used to eliminate the function H from the 

linearized time-dependent Bernoulli equation. The result is 

- (Rtf) ) + RQ(s)(j) = 0)^(l+|Bj)R(i) , (4.4) 

-ns s n 

where 

Q(s) = BR^ + (m/R)^ - + (Z^/R)^] . (4.5) 

s s s 

The solution of differential equation (4.3) gives the boundary 
condition on ^ along the equilibrium surface. The boundary condition 
is specified by the contact angle conditions at s = 0 and s = S , 
which are assumed to be time independent. The perturbed surface and 
the equilibrium surface must have the same contact angle with the 
cylinder wall and the bottom. 

These conditions relate H , ^-3 and X at s = 0 and 

S 3 ^ 

s = S . In terms of the function (j) these conditions are 

sin Y - <}>„ cos y - x^ cos x) = 0 at s = 0 (4.6) 

Xlo XL 0 ^ 

A sin Y + d) Tli cos y = 0 at s = S . (4.7) 

^ns ^n^s 

Equations (4.1), (4.2), and (4.4) - (4.7) determine the eigenfunctions 
(|) and the corresponding eigenvalues u. 

It will be convenient to let 00 , 01, 02 , . . . . denote normal 

modes having cos(m0) dependence with the values m= 0, 1, 2,..., 
respectively. We also shall let RO, Rl, R2, ... denote normal modes 
having 0, 1, 2, ... radial nodes in the interval 0 < s < S (not 
counting the nodes, if any, at the endpoints of the interval). 
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5. Discrete Representation and the Solution of Laplace’s Equation 

The functions (|) and on the boundary of the vertical cross 

section of the liquid will be represented by their values at H + M 
points 

^ = ((f)(Sj^),(f»(s2),. 

and similarly for . These points are shown in Figure 3, The first 
N of these points will be along the meridian of the free surface in the 
cross-sectional plane. The remaining M will be on the meridians of the 
cylinder wall and the bottom in the cross-sectional plane. None of the 
Sj are corner points of the boundary. 

We shall partition the vectors # and into two parts : 

includes values of <|) at points on the free surface, and $2 includes 

those on the cylinder wall and the bottom 

$ = $2> 

~ ( 4 * (® 2 ^ » • * • > 4 * ) 

The boundary condition on (J) along the cylinder wall and the bottom. 

Equation (4.2), becomes 


«2n - “ • 


(5.1) 


Because (j) satisfies Laplace’s equation. Green’s formula yields 

an integral equation that relates ({> and on the boundary of the 

vertical cross section of the liquid. This can be approximated by a , 

matrix equation of the form: 

Wf = Cf . 
n 

The calculation of the matrices W and C is -described in [S]. 


(5.2) 
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r 

Figure 3. The discrete set of points on the boundary. 
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6. Discrete Representation of the Boundary Condition on the Free 
Surface 

Usiug Equation (4.4) and the contact angle conditions (4.6) and 
(4.7), we will derive a, discrete set of equations relating f|) and 
(|!^ at the points . . .3^^, along the meridian of the free surface 

in the cross-sectional plane. Let this meridian be divided into N 
intervals. The interval has 


t. < s < t. , 
3 1+1 


where t, = 0 and t„,, = S . Let s. he the midpoint of the j 
1 iH-1 1 


.th 


interval. 


,th . 


We integrate Equation (4.4) over the j interval 


'^j+1 


R(t.)4> (t.) - R(t.,T)c|) (t.,,) 1- 

j ^ns j 3+1 ^ns 3+I 


Q(s)R(s)<|)^(s) ds 


^^3+1 


0) (1+1b|) 


R(s)^(s) ds 


j 


( 6 . 1 ) 


The integrals are approximated by 


'b+i 


^ (s.) 

a. 1 


Q(s) R(s) ds 


t, 

3 


and 


^(Sj) 


'i+1 


R(s) ds 


t. 

1 
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If iS’ not an endpoint of the meridian, we approximate by 


^ (s ,) - (}) (s . ) 

i Ctj «. la— I !a_jzi- 


ns 3 


s . ~ s . - 
3 J-1 


Substituting these approximations into Equation (6.1) gives 


where 


^^3+1 


A.. - (l+|Bi) 


R(s) ds 


t. 

3 


^^3+1 


P. = 
3 


Q(s)R(s) ds 


Jt. 


and 


■^3,3+1 ^^^3+l^^^®3+l"^j^ 


T, . = - T. , . - T. . + P. 

33 3-l>3 J.J+1 3 


For tj^ = 0 , the inner endpoint of the .-.meridian, we approximate 






=l-‘l 


To eliminate the unknown d) (t.,) , we use the contact angle 

n 1 

condition (4.6) 

-I-i . 


where 




These give 

■ *ns<V ' h • 



127 


where 

/[sin Y + (s^-tj^)L^ ] . 


Substituting this approximation into Equation (6.1) gives 


(6.3) 


where 


"^ll “ ^12 ■ 

For “ S , the outer endpoint of the meridian, we approxi- 


mate 


(j) (t„ . T ) 

^ns N+1 


<I> (t„ - ) - c() (s ) 
^n N+1 ^n n 

^N+1 “ ®N 


To eliminate the unknowi contact angle condition 


(4.7) 


*?> T = “ 1-0 <l> (t^.n) J 

^ns N+1 ' 2 ^n N+1 


where 


h - i’s<W ^ 

These give 

“^ns^^N+l^ “ ^1 *^n^®N^ ’ 


where 

Kj - 1^ / [sin Y + (t^i - 

Substituting this approximation into Equation (6.1) gives 
where 


"^NN .^N-1,N 


^2^^^N+1^ 


+ P 


N 


(6.4) 



128 


Equations (6.2) - (6.4) can be written in matrix form as 

Tf- = 03^' A (6.5) 

in 1 

This is the boundary condition on along the meridian on the free 
surface. A is diagonal, and the diagonal elements are positive. T is 
tridiagonal and symmetric, and the off-diagonal elements are negative. 

The set of Equations (5.1), (5.3), and' (6.5) is the discrete version 
of the eigenvalue problem for the small- amplitude, periodic sloshing 
modes of a liquid in a vertical, rotationally symmetric cylinder. 
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7, numerical Solution, of tlie Discretized Eigenvalu'e Problem 

We write the matrices W and C of equation (5.3) in block form; 
.W^^, ,Wj^ 2 » ^2V' ^22’ similarly for C. Subscript 1 denotes the rows 
and colximns corresponding to the N points along the free surface, and 
subscript 2 denotes those corresponding to the M points along the- 
cylinder wall and- bottom. Since $ 2 ^ zero, Equation (5,3) can be 
written 


W 0+W $=G $ 1 

'^11 ^1 *^12 2 11 In ’ 


TJ <1 + W # = $ 

r «22 ^21 ^In 


(7,1) 


The matrix A is diagonal, so Equation (6,5) is easy to solve for , 
which we can eliminate from Equations (7.1), 

Define 

¥ = ¥ A T 

11 11 ’ 


F = W A ^ T 
21 21 ' 


Then Equations (7.1) give 


(7.2) 


^11 ^In “ ^^ll®ln “ ^2^2^ ’ 


^21 ^in = ^ (^2l"ln ‘ ^22^2^ 


(7.3) 


Equations (7.3) can be written as single matrix equation for the 
eigenvector 


11 


'21 



to 



ll 


“^22/ \^2 


(7.4) 


Equation 7.4 could be solved for the eigenvalues w j however, M of 

. 2 

the eigenvectors have the eigenvalue a) = 0 . A linearly independent 
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set of these eigenvectors is 
* 1^=0 

where e^ is the vector with a one in the position and zeros else- 
where. These eigenvectors correspond to no motion of the free surface, 
since is zero. A computer program that calculates all the eigen- 

values of a matrix, such as the IMSL routine EIGZF, will waste some time 
computing these unwanted eigenvalues. 

We can avoid calculating the zero eigenvalues by eliminating 
from the pair of Equations (7.3). Define 


D = Cii 

-1 

-"11 ”22 

*^21 ’ 


E = Fii 

-“12 "S 

*■21 ' Wll - ”12 ”22 

(7.5) 


Then Equations (7.3) combine to give 


Ei>. = 0) 

In In 


(7.6) 


Equation (7.6) can be solved for its eigenvalue by the IMSL routine 
EIGZF, which uses a QZ algorithm to reduce E to upper Hessenberg form 
and D to upper triangular form. 

The solution of Equation (7.6) is performed by the program SLOSH. 
The input to SLOSH is the set of points describing the equilibrium 
free surface calculated by CAPIL and parameters that define the 
cylinder wall and .spheroidal bottom. SLOSH then calculates the 
matrices A, T, W, and C; uses the IMSL routine LINVIF to calculate 
W 22 inverse; calculates the matrices D and E; and uses EIGZF to 
calculate the eigenvalues. This method of solving the eigenvalue 
equation is not the most computationally efficient, but by using 
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the existing and reliable IMSL routines it requires the least amount 
of programming effort. 

For comparison, the routine EIGZF was used to solve both Equation 
(7.6) and Equation (7.4) for a few cases,. The numerical values of 
corresponding eigenvalues for these two methods were identical to the 
four figures that were printed out in each case. 
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8 . Small" Amplitude Periodic Sloshing Modes of the Liquid between Two 
Concentric Right Circular Cylinders 
In this section we solve the eigenvalue problem for the small- 
amplitude, periodic sloshing modes of the liquid contained between 
two concentric, vertically oriented, right circular cylinders of radii 
Tq and r^^ . A cross' section is shown in Figure 4. The equilibrium 
surface is a horizontal plane when the contact angle is 90*^. The nor- 
mal mode problem for this case has an analytic solution. We can use 
this solution to test the accuracy of the program SLOSH. 

Let V be the rectangular domain r^ < r < r^ and 0 < z < z^ . 
Laplace’s Equation for ({) in the domain V is 


(j) +-(j)+(|) -2L-(}) = 0 

Yrr r ^r ^zz 2 ^ 

r 


The boundary conditions are 


( 8 . 1 ) 


d) = 0 
^r 


at r = r^ and r^ 


d) = 0 
^r 


at z = 0 . 

Equation (4.4) for (j) on the free surface, z 




( 8 . 2 ) 


where 

Q(r) = B + (m/r)^ 

The contact angle conditions. Equations (4.6) and (4.7), become 

4>^z ^0 

We solve these equations by separation of variables. Let 
4>(r,z) = X(r) U(z) 
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Figure 4. Cross section of two concentric cylinders or radii ro 


and ri , respectively. 
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then (8.1) gives the pair of equations . 


U" = U (8.3) 

X" + i X' + (k^ - m^/r^)X = 0 , (8.4) 

with boundary conditions 

U’ = 0 at z = 0 , (8.5) 

X' = 0 at r = Tq and r^^ (8.6) 


The contact angle conditions will be automatically satisfied if Equation 
(8.6) is satisfied. 

Equation (8.2) gives an -equation for the eigenvalue 

Ofeg) 

The solution of (8.3) and (8.5) is 
U = cosh (kz) 

Thus the eigenvalue is 

2 k(k^ + B) 

= (1 +-1b|) tanh(kZp) . (8.7) 

The solution of (8.4) is 

X(r) = c J (kr) + d Y (kr) , 
ra m 

where J and Y are Bessel functions of the first and second kind of 
m m 

order m. Equation (8.6) requires 

c J ’ (kr„) + d Y ' (kr„) = 0 , 
m 0 m 0 

c j ' (kr., ) + d Y ' (kr ) = 0 , 

m J. ml 

These will have a nontrivial solution for c and d if 



( 8 . 8 ) 



Equation (8.8) gives the values of k for the normal modes. The first 


few values for the case r^ = 0.5 and r^ - 1.0 are listed in Table 
1 to the accuracy indicated. 


Table 1. k values for m = 0, 1, and 2. 


m = 0 

‘ m = 1 

m = 2 

0.0 

1.3547 

2.6812 

6.3932 

. 6.5649 

7.0626 

12.6247 

12.7064 

12.9494 

18.8889 

18,9427 

19.i032 

25.1624 

25.2045 

25.3224 


The solution m = 0, k = O'. 0 corresponds to no movement of the 
equilibriimj surface or of the liquid. 
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9. Comparison of the Analytically and Numerically Calculated Solutions 
for the Normal Modes of the Liquid between Two Cylinders 

Figure 4 shows the cross section of a liquid contained between two 
concentric right circular cylinders, oriented vertically. Each of 
the four sides of the cross section (the free surface, the bottom, and 
the two cylinder walls) was divided into n intervals of equal length. 
The velocity potential the perimeter of the cross section was 

represented by its values at the midpoints of these 4n intervals. 
These 4n values of (|> are related by Equations (5.1), (5.3), and 
(6,5). Nianerical solutions of these equations were computed for the 
case Tq = 0.5, r^^ = 1.0, Zq = 0.9 , contact angle 90°, and Bond 
number 0 using the program SLOSH. 

Numerically calculated squares of the frequencies for the modes 
01RO, 01R1, and 01R2 using n = 5, 10, and 20 points are shorn in Table 
2. The corresponding analytic values- for the squares of the frequen- 
cies, calculated from Equation (8.7) and the k values of Table 1, are 
also shown. 


Table 2. Squares of frequencies for various normal modes 



01RO 

eiRi 

01R2 

5 points 

2.157 

292.1 

1965. 

10 points 

2.127 

284.8 

2080. 

20 points 

2.118 

280.7 

2072. 

analytic 

2.087 

282.9 

2052. 
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The relative errors of the numerically calculated squares of the fre- 
quencies [(jd (n points) - (analytic)]/ tii (analytic), are shown in 
Table 3. 


Table 3. 

Relative errors 

of the squares 

of frequencies. 


01RO 

01R1 

01R2 

5 points 

0.034 

0.033 

-0.042 

10' points 

0.019 

0.007 

0.014 

20 points 

0.015 

-0.008 

0.010 


Note that the relative errors of the frequencies are approximately half 
these values. The error decreases substantially between n = 5 and 
n = 10, but less so between n = 10 and n = 20. Even the errors for 
n = 5 are quite small, considering that only five radial modes can be 
represented by a 5-point approximation to the meridian on the free 
surface. 

Ntimerically calculated squares of the frequencies for the modes 
0ORO, 0OR1, 0OR2, 02RO, 82R1, and 02R2 using n = 10 points are shown 
in Table 4. Corresponding analytic values are shown also. 


Table 4. Squares of frequencies for various normal modes. 



0ORO 

0OR1 

0OR2 

10 points 

0.7 10“^^ 

265.4 

2038. 

analytic 

0.0 

261.3 

2012. 


82R0 

02R1 

02R2 

10 points 

19.99 

346.6 

2201. 

analytic 

18.97 

352.3 

2171. 
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The relative errors of these squares of frequencies are shown in Table 

5. 


Table 5. 

Relative Errors 

of the squares 

of frequencies . 


RO 

R1 

R2 

00 modes 

— 

0.016 

0.013 

02 modes 

0.054 

-0.016 

0.014 


The relative errors of the squares of frequencies for n = 10 points, 
as shoxm in Tables 3 and 5, are typically from 0.01 to 0.02. These 
results show the program SLOSH calculates with satisfactory accuracy 
for our purposes the frequencies of the normal modes of a liquid con- 
tained between two concentric right circular cylinders. 
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10. Equilibrium Free Surfaces of a Liquid In a Vertical Right Circular 
Cylinder with a Concave Spheroidal Bottom 

With a given volume of liquid in the cylindrical container we 
associate a dimensionless fill height defined as follows: let the 

given volume V equal the volume bounded by the container wall and 
bottom and the horizontal plane z = z^ . Then the fill height h^ 
is z divided by the container radius a. 



The axial ratio of the spheroidal bottom is b/a = 0.724 . 

Equilibrium free surfaces, approximated by 21 points on the meri- 
dian, were calculated by the program CAPIL for contact angle y = 0° 
and for the fill heights: 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 

0.60, and 0.70. For each fill height equilibrium surfaces were 
calculated for a sequence of increasingly negative Bond numbers. The 
first surface for each -fill height was calculated for Bond number 
B = 0. The initial approximation to this surface was generated by the 
subroutine CYLCUR. The equilibrium surface for each Bond number was 
used as the initial approximation to the surface for the next Bond 
number in that sequence. 

The equilibrium surfaces that we .have calculated are members of a 
family with two parameters B and h^ . Let B^^ denote the critical 
value of the Bond number for the nonexistence of equilibrium surfaces 
of this family for a given fill height. Let denote the criti- 

cal value of the Bond number for the stability of equilibrium surfaces 
of this family for a given fill height. Stable equilibrium surfaces 
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exist for B . B , unstable equilibrium surfaces exist for 
crxt 

B '^B<B , <0ifB and no equilibrium surfaces of 

eq crit eq crxt ; 

this family exist for B < B . (Other equilibrium surfaces might 

eq 

exist, such as multiple-valued surfaces or surfaces with shapes very 
different' from those of this family.) Concus and Karasalo showed that 
unstable equilibrium surfaces exist for B infinitesimally lower than 

a. 

B and h < h*' = 0.503 , but that no equilibrium surfaces of this 

crxt V V — 

family exist for B < B . and h ^ h [l] . Their result may be 

crxt V V 

restated as B < B . for h < h , but B = B .. for h ^ K 

eq crxt v v eq crxt v v 

Our calculations agree with their result and provide an estimate of 



For each fill height we found a Bond number » depending on 

h^ , for which the iteration for the equilibrium surface diverged. 

The iteration using was approached' by a sequence of calculations 

using small decreases in B . Let B denote the Bond number imme- 

diately preceding that sequence, ®conv ® 

B ^ ®eonv equilibrium surface changed only slowly with B. The 

equilibrium surface for each value of B was an excellent approximation 

to that for the next value of B in the sequence. This indicates that 

the divergence for the case was caused not by the initial 

approximation but by the nonexistence of an equilibrium surface for 

this family. Thus approximation to B^^ . Table. 6 shows 

B and B,. as a function of h . It also shows B .. calcu- 

conv dxv V crxt 

lated to four decimal places by Concus and Karasalo [l]. 



141 


Table 6. 

B , B , . , and 

conv dxv 

B for various 

crxt 

fill heights. 

h 

V 

~B. 

conv 

dxv 

-B 

crxt 

0.20 

1310. 

1320. 

480.4283 

0.25 

488. 

492. 

238.6539 

0.30 

216. 

218. 

•132.9638 

0.35 

107. 

108. 

79.6741 

0.40 

58.0 

58.2 

49.9096 

0.45 

33.4 

33.6 

31.9190 

0.50 

20.2759 

20.2760 

20.2759 

0.60 

8.42 

8.43 

8.4411 

0.70 

3.98 

3.99 

4.0020 


The data of Table 6 are shown in Figure 5. The solid* line is the 
graph of B .. and the dashed line is that of B, . . These lines 
divide the Bond-number, fill-height parameter space into three regions: 
one for which stable equilibritmi surfaces exist, one for which unstable 
equilibrium surfaces exist, and one for which no equilibrium surfaces 
exist. (Growth rates for perturbations of the unstable equilibrium 
surfaces were calculated by the program SLOSH for various values of B 
and h^ . These will be discussed in the next section of this report.) 

Figure 5 also shows data points from stability experiments carried 
out at the NASA Lewis Zero Gravity Facility for the container shoxm in 
Figure 1 [6]. The experiments used three containers with radii 7 cm , 
5.5 cm' , and 2 cm', respectively. In the experiment the container had 
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Fieure 5. B . and B,. as functions of fill height compared 
® crxt div 


with experimental points. 
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approximately 2.5 sec .of free fall followed by approximately 2.5 sec 
of negative low-g fall. During the first 2.5 sec the liquid surface 
adjusts from one g to zero g . During the next 2.5 sec instabil- 
ities may be observed if they grow sufficiently rapidly. Solid data 
points correspond to experimental parameter values for which the 
surface was observed to be unstable. Open data points correspond to 
parameter values for which the surface did not develop a noticeable 
instability within the 2.5-sec time interval. The experimental data 
and the n-umerically calculated curves agree quite well. All the 
experiments in which the surface was observed to be unstable have Bond 

numbers B' < B , <- 0 . 

crit 

®div approximation to B^^ , The accuracy of this approxi- 

mation can be investigated by considering the cases h^ = 0.60 and 
0.70 . For - 0.60 CAPIL diverges for some Bond number in the range 

(-8.42, -8.43) and- B = -8.4411 . The relative error in this case is 

eq 

less than 0.003. For h = 0.70 CAPIL diverges for some Bond number 

in the range (-3.98, -3.99) and B^^ = -4.0020. The relative error in 

this latter case is less than 0.005. 

From this we infer that the correct value of B for h =0.50 

eq V 

is slightly less than the value found here for the 21 point surfaces, 

** M 

and that there is a small range of Bond numbers between B . and 

crit 

B^g for this- value of h^ . This is supported by a calculation of 
the frequencies of individual normal modes, which is discussed in the 
•next two sections. Based on an approximate calculation of the fre- 
quencies, the RO01 mode becomes marginally stable at B = -20.243 , 
while all the other modes approach instability as ‘B approaches -20.276. 
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, 11. Frequencies of the Normal Modes of a Liquid in a Vertical Ri^ht 
Circular Cylinder with a Concave Spheroidal Bottom 
The frequencies of the small— amplitude periodic sloshing modes of 
a liquid in a vertical right circular cylinder with a concave spheroidal 
bottom were calculated by the program SLOSH for contact angle Y “ 0° • 
The axial ratio of the spheroidal bottom is b/a = 0.724 . The equi- 
librium free surfaces were approximated by 21 points of the meridian, 
as described in Section 10. These 21 points were the endpoints and 
midpoints of 10 intervals on the meridian. The velocity potential 

for perturbations to these surfaces was represented by its value at 
the 10 midpoints of these intervals, by its value at 10 points on the 
meridian on the cylinder wall, and at 10 points on the meridian on 
the bottom. We shall refer to this as the 10-point approximation to 
(|) . Surfaces corresponding to numerous values of h^ and B were 
used. 

A few surfaces approximated by 41 points on the meridian were 
used to check the accuracy of the frequencies calculated using the 21 
point surfaces. For these cases (j) was represented by its value at 
20 points each on the meridians on the free surface, the cylinder 
wall, and the bottom. We shall refer to this as the 20-point approxi- 
mation to (j) . 

The squares of frequencies for various normal modes and for vari- 
ous values of h^ and B calculated by SLOSH using the 10-point 

approximation are shown in Tables A 1 through A 11. Typically the 
2 

values of U) in these tables have a relative error of 1-2% for 
2 

values of w that are not too small and for Bond numbers that are not 
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too near B,. . This will be discussed in more detail in Section 12, 

dxv 

These squares of frequencies are plotted as functions of B in 

Figures 6-12 for h^ = 0.20, 0.30, ...0.70 . Note that in Figures 

2 

10-12 ( h = 0.50 - 0.70 ) a different scale for O) is used for each 

V 

mode that is plotted-. The purpose is to show that all these modes have 

2 

a similar dependence of to on B . However, in Figures 6-9 ( h^ = 

0.20 - 0.40 ) all the RO modes ( RO01, R092, RO03, ... ) that are 

2 

plotted in a given figure use the same scale for (O . This is to show 

for each value of B which mode is most negative. 

We shall first describe the general features *of these figures, 

and then consider numerical details for particular cases and discuss 

the accuracy of the calculations. 

2 

Figure 6 shows graphs of to- (B) for the modes R061, RO02, RO03, 

RO04, RO06, and R10O for h^ = 0.-20 , The mode R10O is plotted with 

a scale 1000 times that of the other modes. (An accuracy check shows 

the values for the R10O mode are about 10% too large. However, we 

include it in Figure 6 for a rough comparison with the RO modes . ) Note 

that the Bond numbers for which the various RO modes become marginally 

2 

stable [for which 0 ) (B) = O] . lie in a small range. 

Figure 7 shows -this more clearly. The order in which the modes 

become unstable is as one would expect: First RO01 , then RO02 , 

RO03 , RO04 , and RO06 . (RO05 was not calculated for h^ = 0.20 .) 

Note also that because the higher '0 modes have steeper slopes 
2 

(do) /dB ), each mode, in turn, becomes the dominant unstable mode 

2 • 

(most negative value of to ) for a short range of Bond numbers. 
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Bond number 


2 

Figure 6. 0 ) (B) of various modes for h = 0.20 . 

V 
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This same pattern for the RO modes is shown in Figure 8 for 

2 

h^ = 0.30 . A new feature appears in this figure. It is that w 
for the R10O mode passes through a point of inflection and begins to 

curve downward,. The rate at which it approaches zero, the magnitude of 

2 2 
do) /d(B) , increases as B approaches . The functions co (B) 

for the other R1 modes, R101 to R106 , have similar shapes and 

differ from that of R10O by only a few percent, as shown in Table 

2 

A 3. The functions to (B) for the modes R20O to R206 have shapes 
similar to those for the R1 modes but magnitudes about five times 
larger. All the R1 and R2 modes curve downward as B approaches 


®div 

2 

Figure 9 shows that for h^ =0.40 OJ becomes negative between 
B . and B only for four modes: R061 , RO02 , RO03 , and 

RO04 . Furthermore, only the first three modes become dominant insta- 
bilities in this range. The rate at which the R10O mode approaches 


2 

zero, du /dB| , becomes very great as B approaches B 


div 


Note 


also that u'^(B) for each of the RO modes passes through an inflec- 
tion point and curves downward as B approaches • 

Figures 10, 11, and 12 are for the cases h^ = 0.50 , '0.60' , and 

2 

0.70 , respectively.. In each case the functions u (B) for the vari- 
ous modes have similar shapes. They all curve downward for B near 
B^rit » rate at which they approach zero becomes very great 

as B approaches . All modes apparently go to zero at or near 

2 

^crit ’ "Bhis behavior of w (B) is consistent with the nonexistence 

A 

of an equilibrium free surface -nearby the critical one for h_^ ^ 
and 'B < < 0 . It is in sharp contrast to the behavior seen in 







XBL78M00 
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ft)2(B) for RO mode 
co%B)/5 for RO 62 mode 
co 2(B)/15 for RO 63 mode 
cy2(B)/50 for R1 60 mode 
fi;2(B)/100 for RO (96 mode 

0 


Bond number 

2 

Figure 12. (o (B) of various modes for li = 0.70 . 
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the' cases h = 0.20-0.40 , in which only a few RO modes were un- 

stable for a range of Bond numbers beyond . 

For the case h^ = 0.50 the RO01 mode becomes marginally stable 

at a slightly higher Bond number than the other modes. In the 10- 

point approximation it becomes marginally stable at B = -20.243 , 

while all the other modes approach instability as B approaches 

-20.276 . It appears that for this case there exists a very small 

range of Bond numbers between B _ and B 

crxt eq 
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12. Frequencies of the Normal Modes Continued — Accuracy 

2 

The bulk of our data are values, of co calculated by the 10-point 

approximation. Tables A 1 - A 11. Throughout this section we shall 

2 

investigate the accuracy of these data. The few values of o) calcu- 
lated by the 20-point approximation are used solely to estimate the 

accuracy of these -data and how they can be Improved. We shall refer 

2 2 

to these values as ( 0 ^^ and 0 ) 2 q > respectively. 

2 

We are most interested in the accuracy of 0 )^^ for the growing 

2 

RO modes. These negative values of to are small numbers, so a 

small (absolute) error in them can be significant. We shall show that 
2 2 

the error in {Oj^q(B) is approximately -(dto /dB)AB* , where AB* is 
a function of h but hot of the mode number, that is, 

(0^(B) (dto^/dB)AB* , (12.1) 

or, equivalently, 

tO^(B-AB*) « • C12.2) 

2 

Thus, a value of (o^q(B) from Tables A 1, A 2, or A 4 actually corres- 
ponds to the Bond number B-AB* . 

Five comparisons of calculated quantities support this description 

2 

of the approximate dependence of the error in oj^q on the parameters 
h^ , B , and mode number. The first comparison is between the values 
of calculated by the 10-point approximation, which we shall 

denote by those calculated to four decimal places from 

a static analysis [l]. AB* is defined as the difference in these 


values . 
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®crit,10 " \rit ' 

These quantities are shown in Table 7. 

of B _ . 
crxt 


AB* is ap.p,roximately 1 - 2% 


Table 7. B B three fill heights 

crit critjlO 


h 

V 

-B 

crit 

®crit,10 

Ab* 

0.20 

480.43 

468.2 

12.2 

0.30 

132.96 

130.73 

2.23 

0.40 

49.91 

49.35 

0.56 


2 

The second comparison is between the values of o) (^^rit^ 

RO01 mode calculated by the 10-point and 20-point approximations. 

2 

Since the correct value of (i) is zero in this case, these values 

2 

are errors. They show that the error in the calculated values of a) 

2 

depends as 1/N on the number of points used to approximate <j> . 

2 

These values are shoxm in Table 8. Note that W 2 Q is approximately 
1/4 of • 


? 2 

Table 8. and a)„„ at B .. for the RO01 mode. 

10 20 crit 


h 

V 

2 

“10 

2 

■““20 

0.20 

0.00180 

0.00045- 

0.30 

0.00183 

0.00046 

0.40 

0.00209 

0.00053 
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2 2 

The third comparison is between and for various RO 

2 

modes, fill heights, and Bond numbers. Define Aw as 




2 ■ ' 2 
Since the error in is approximately 1/4 of the error in , 

2 2 
it follows that Aco is approximately -3/4 of the error in 0 )^^ . 

2 

The values of Aw are shown in Table 9. They vary greatly with 
mode number - 


Table 9. 

2 

Ato for various 

RO modes , 

fill heights, and 

Bond numbers. 

h 

V. 

-B 

R091 - 

R092 

R083 

0.30 

130. 

0.00141 

0.0055 

0.0120 


132.96 

0.00137 




140. 

0.00127 

0,0050 

0.0108 

0.40 

45. 

0.00193 

0.0070 



49.91 

0.00156 




50. 


0.0059 



Define AE as 
AB = Am^(dB/dm^) 

2 

The values of dci) /dB can be calculated approximately by central 

differences of the data in Tables A 2 and A 4. The resulting values of 

AB are shown in Table 10. Mote that for a given fill height, while 
2 

Aw varies greatly with mode number, AB does not. Our fourth com- 
parison is between the values of AB and AB* . Table 10 shows that 
the values of AB are approximately 3/4 of the corresponding values 
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2 2 
of AB* . Since Au is approximately -3/4 of the error in , 

2 2 

this implies that the error in is approximately -(dw /dB)AB* 

for some range of Bond numbers containing • 


Table 10. AB for various RO modes, fill heights, and Bond numbers. 



For a given fill heights, let Bg^ denote the Bond number for 

which the R09n mode becomes neutrally stable, that is, for which 

w (Bg^) = 0 . Table 11 shows the values of the Bond numbers for the 

neutral stability of various modes as calculated by the 10-point 

approximation to 4> , which we shall denote by B„ • As B 

bn,±U 

increases the various modes become unstable in order of increasing 

9 mode number, so Bg^ is . We have already compared the 

accurate values of B . ^ with the corresponding values of B , » 

crxt ® crxtjlO 

91,10 ■ 


that is, with B 
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Table 11. -Bp, for various modes and fill heights. 


h 

V 

Boei 

RO02 

R093 

RO04 

RO05 

RO06 

0.20 

468.25 

471.42 

475.42 

480.89 


495.85 

0.30 

130.73 

133.22 

137.11 

142.49 

149.28 

157.60 

0.40 

49.35 

51.28 

54.19 

57.19 




Our last comparison is between the values of and 

'in • Define AB- as 
9n,20 9n 

"" ®0n,lO ” ^6n',20 * 

Table .12 shows values of ABg^ for various modes and fill heights. It 

shows that the values of AB. are approximately 3/4 of the correspond- 

Un 

ing values of AB* . For a given fill height approximately 

2 

the same for each mode. This implies that the error in is 

2 

approximately -(dw /dB)AB* for some range of Bond numbers containing 

these B . 
on 


Table 12. 

ABo for 

0n 

various modes 

and fill heights. 

h 

V 

RO01 

■RO02 

RO03 

0.30 

1.68 

1.69 

1.59 

0.40 

0.42 

0.42 



When the values of Bg^ given in Table 11 are adjusted by adding 


AB* , the values shown in Table 13 are obtained. 
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Table 13: 

Adjusted values of 

-®6n 

various 

modes and 

fill heights. 

h 

V 

.RO.ai. . 

RO02. 

RO03 

■ RO04 

. RO0-5 

RO06 

0.20 

480.43 

483.6 

487.6 

493.1 


508.0 

0.30 

132.96 

13 5. '5 

139.3 

144.7 

151.5 

159.8 

0.40 

49.91 

51.8 

54.8 • 

57.8 




We consider finally the R1 modes. Table 14 shows the relative 
2 

difference in 0 ) calculated with the 10-point and 20-point approxi- 

2 2 2 2 

mations, that is, C^20 ~ '^'lO^^^lO * value of 0) differed by 10% 

for the R101 mode with h = 0.20 . It differed by 2.6 - 2.7% for the 

V 

R101 and R102 modes with h =0.40 for Bond number 55, which is 

V 

2 

• near • The value of a) differed by only 1-2% for all the 

R1 modes with h =0.30 or 0.40 and Bond numbers not near B,, 

V div 


2 2 

Table 14. Range of Ato /w for various modes and fill heights. 


h 

V 

range of -B 

R101 

R102 

R103 

0.20 

450.-500. 

0.100 



0.30 

130.-140. 

0.015-0.016 

0.015-0.016 

0.015-0.016 

0.40 

45.- 50. 

0.015-0.019 

0.014-0.019 


0.40 

55. 

0.026 

0.027 
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13. Growth Rates and Accuracy for Fill Height = 0.30 

2 

Figure 8 and Table A 2 show oj (B) of the various RO inodes for 
h^ = 0.30 . They show which is the maximally unstable mode for each 
value of B . The information for .the maximally unstable mode is 
.displayed in Table 15. The values of B listed in this table have 
been adjusted by AB* . 

■For example, the RO01 mode is the maximally unstable one for • 

-134.04 < B < -132.96 , the RO02 mode for -140.27 < B < -134.04 , 

2 2 
etc. The value of (du /dB)AB* is an estimate of the accuracy of u 

2 

before adjustment by AB* . We feel the error remaining in (O after 

2 

this adjustment is less than (dm /dB)AB* . In particular, we -feel 

2 

the errors remaining -in w for the RO01 , RO02. , and RO03 modes 

are 1/10 to 1/4 of (doj^/dB)AB* . 


2 

Table 15. w (B) of the maximally unstable mode for h^ >= 0.30. 


-B 

2 

0) 

maximally 

unstable 

(do)^/dB)AB* 

132.96 

0.0 

Roei 

0.0017 

134.04 

0.0025 

RO02 

0.007 

140.27 

0..0218 

RO03 

0.014 

149.99 

0.0808 

RO04 

0.021 

163.48 

0.204 

RO05 

0.028 

179.76 

0.407 

RO06 

0.032 

202.23 

0.701 
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The dimensionless growth rate F of the maximally unstable mode 
is shown in Table 16. The corresponding growth, period in seconds is 

■r“^.(t/t) = r“^[pa^/a(l+lB|)]^ . 

This is calculated for a cylinder of radius 7 cm for the three 

liquids ethanol, freon, and FC78. The values of p/a used for these 

2 3 

were 0.03538, 0.08489, and 0.131 sec /cm , respectively. Ethanol has 
the fastest growth rates and FC78 has the slowest. -'132.96 

for this case. At Bond number B = -150 the growth periods range 
from 1.0 to 1.9 sec. At B = -202 , which is 50% beyond j they 

range from 0.29 to 0.56 sec. It is not likely that growth would be 
observed in these cases in an experiment with a negative-B phase of 
only 2.5 sec, since only 2-8 growth periods would elapse. 


Table 16, Maximal growth rates and growth periods for ~ 0.30 . 


dimensionless 

values 

growth period 

(sec) 

-B 

r 

ethanol 

freon 

FC78 

132.96 

0.0 

CO 

00 

CO 

134.04 

0.050 

5.9 

9.2 

11.4 

140.27 

0.148 

1.99 

3.1 

3,8 

149.99 

0.284 

1.00 

1.54 

1.92 

163.48 

0.451 

i 0.60 

0.93 

1.16 

179.76 

0.638 

0.41 

0.63 

0.78 

202.23 

0.837 

0.29 

0.45 

0.56 


2 2 
The errors in w for the smaller values of to (the RO01 , 

RO02 , and R.O03 modes) have a greater percentage reduction from the 

2 

AB* adjustment than those for larger values of w . However, these 
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-errors were initially larger fractions of their values of u than 

2 

those for the larger values of u . As a result of these two effects, 
2 

the errors in OJ remaining after this adjustment probably lie in the 

2 

range 5 - 20%, the larger values of u ' being more accurate. The 
corresponding errors in the growth rate probably lie in the range 
2 10% . However, these are only the computational errors in T ; 

they represent the accuracy with which the growth rates were calculated 
from the assumed model of the liquid motion. The accuracy with which 
they describe experimentally observed growth rates depends also on 
the accuracy of that model. In this model the fluid motion was assumed 
to be nonviscous and irrotational, and the contact angle was assumed to 
be time independent. These asstanptions could be tested by computing 
the fluid motion with a complete hydrodynamics code that includes all 


the relevant effects. 
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14, Summary 

In this paper we calculate the small-amplitude periodic sloshing 
modes of a liquid in a vertical right circular cylinder with a concave 
spheroidal bottom, for the case in which there is not sufficient liquid 
to cover the bottom entirely. Numerical results are obtained for a con- 
tainer currently used for the storage of liquid fuels in the Centaur 
space vehicles, for which the axial ratio of the bottom is b/a - 
0.724 . 

We follow the derivation in [2] for obtaining the equations of 
motion for the case studied here, but we use a different technique for 
obtaining the numerical solution. The liquid is subject to surface 
and gravitational forces. The equilibrium surface is the solution of 
the. time-independent Bernoulli equation subject to a contact-angle 
condition. 

It is assumed for the dynamical equations that the fluid flow is 

irrotational and incompressible. The fluid velocity is the gradient of 

a potential function that satisfies Laplace’s equation. The velocity 

potential and its gradient on the free surface are related by the 

linearized time-dependent Bernoulli equation and the contact-angle 

condition. The sloshing motion is analyzed in terms of normal modes. 

The discrete form of these equations yields a generalized eigenvalue 
2 

problem for co , the square of the normal-mode frequency. This 
problem was solved numerically using the IMSL routine EIGZF.. 

The accuracy of this numerical procedure was tested by calculating 
the eigenvalues and eigenvectors for the small-amplitude periodic 
sloshing modes of a liquid contained between two concentric vertical 
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circular cylinders for contact angle y “ comparing with the 

2 

known analytic solution for this case. The numerical values of a) 
were correct typically to about 1 or 2%, a satisfactory accuracy for 


our purposes. 

Equilibrium surfaces of a liquid in a vertical circular cylinder 
with a concave spheroidal bottom were calculated for contact angle 
y == 0" , axial ratio of the spheroidal bottom b/a = 0.724 , fill 
heights h^ ranging from 0,20 to 0,70, and' many values of the Bond 
number. These equilibrium surfaces are members of a family with 


parameters B and h 


B _ was defined above as the critical 
crxt 


value of the Bond number for the stability of surfaces of this family 

for a given fill height. B was defined as the critical value of 

the Bond number for the nonexistence of equilibrium surfaces of this 

family. Stable equilibrium surfaces exist for ^ B , unstable 

equilibrium surfaces exist for B <B<B ..,<0 if B ?^B 

eq crit eq crit 

and no equilibrium surfaces exist for B < B 

eq 

For all the values of the fill height that were studied, stable 
equilibriinn surfaces were found for a range of Bond numbers, 

^crit ^ ® ® accuracy of these calculations, we found the 

same value for B as was found in the static analysis of the same 

problem [l3 . 

For fill heights ranging from 0.20 to 0.45, we found unstable 


equilibrium surfaces for a range of Bond numbers B B < B _ , 

° conv crit 

but no equilibrium surfaces of this family were found for 

B < B,. < B . ( B and B,. are approximations to B .) 

div conv conv div eq 
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I'or = 0.50 unstable equilibrium surfaces were found for a very 
small range of Bond numbers. For = 0.60 and 0.70 no equilibrium 


surfaces of this family were found for B < B 


crxt 


To the accuracy 


of these calculations, these results are consistent with [l], which 

fotmd that B = B .. for h > h* == 0.503 , but that B < B , 
eq crit v v eq crit 

for h < h* . 

V V 

The qualitative nature of the stability of the individual normal 
modes differs for the two cases h < h* and h > h* . For fill 

V V V V 

heights h^ = 0.20, 0.30, and 0.40 , the normal modes RO01, RO02, 

RO03, ... become marginally stable at a sequence of Bond numbers, 

... B „ < B „ < B , = B < 0 . Each RO mode is the fastest 

growing mode for a small range of Bond numbers. For fill heights 

h^ = 0.60 and 0.70 all the modes that were studied approach 

instability as the Bond number approaches each mode the 

2 2 

function o) (B) curves toward the u = 0 axis, approaching zero 


with increasing rapidity as B approaches B 


crit 


For h = 0.50 , 

V 


which is near the critical fill height h* , the RO01 mode becomes 

marginally stable at a slightly higher Bond number than the other 

modes. The instability of all modes for h > h* and B < B 

is consistent with the nonexistence of equilibrium surfaces nearby 

the critical one for this range of parameters. 

2 

Most of the calciilations of to were made by approximating the 

velocity potential on a meridian along the free surface by its value 

at 10 points. It was possible to correct partially these calculated 
2 

values of to by applying an adjustment based on a study of the 


errors . 
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Growth rates of the maximally -unstable mode were calculated for 

2 

the case h^ = 0.30 using the adjusted values of w , Each of the 
modes RO01, R,O02, ... RQ06 , in succession, was the maximally unstable 
one for a small' range of Bond numbers. The corresponding growth periods 
in seconds were calculated for a cylinder of radius 7 cm for the 
three liquids ethanol, freon, and FC78. Ethanol has the fastest 
growth rates and RG78 has the slowest. -132.96 for this 

case. At Bond number B = -150 the growth periods range from 1.0 to 
1.9 sec. At B = -202 , which is 50% beyond > they range from 

0,29 to 0.56 sec. It is not likely that growth would be observed in 
these cases in. an experiment with a negative-B phase of only 2.5 see, 
since only 2 to 8 growth periods would elapse. 
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LIST OF SYMBOLS 


a 


Radius of cylindrical container and horizontal 
semiaxis of spheroidal bottom. 


A 


Diagonal matrix in the discretized time-dependent 
Bernoulli equation. 


b 

B 

B ... 

crit 


Vertical semiaxis of spheroidal bottom. 

2 • 

Bond number = Ka 

Critical Bond number for stability of equilibrium 
surfaces. 



B , B.. 
cony div 

®0m 


Critical Bond number for the nonexistence of 
equilibrium surfaces of the family considered 
in this report. 

Approximations to B^^ 

Critical Bond number for stability of the RO0m 
mode. 


c,d Constants iii a linear combination of Bessel 

functions . 

C Matrix in the discrete solution of the Laplace 

equation. 



171 


D, E, F 

V 

e. 

J 


g 


h 

V 

h* 

V 


1 


H 

H 


0 


H(s,9,t) 

H(s) 

J 

m 


Matrices . 

A rectangular domain. 

Vector with a one in the j position and zeros 
elsewhere. 

Acceleration due to gravity, considered positive 
when directed vertically downward. 

Dimensionless fill height. 

Critical h^ for existence of unstable equilibrium 

surfaces of the family considered in this report. 

These exist for h < h* and B < B < B 

V . - V eq crxt 

Mean curvature at a point on the free surface, 
considered negative when the surface is concave 
upward . 

Scaled mean curvature = Ha . 

A constant = (p^-pQ)a/2a , Interpreted as the 
extrapolated value of H at the height z = 0 . 

Displacement, n of the free surface. ■ 

A factor in the normal mode expression of the' 
displacement of the free surface. 


Bessel function of the first kind of order m . 
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k. 

K, L 

m 

M 


n 


subscript n 
N 


PO 

P. 

3 

Q(s), q(r) 


r 


r 


subscript r 


The argument of the Bessel function is kr . 

Terms representing the contact-angle conditions. 

Number of angular nodes in the normal mode. 

The meridians of the cylinder wall and spheroidal 
bottom in the cross-sectional plane are divided 
into M intervals . 

Each of the meridians of the free surface, cylinder 
walls, and flat bottom in the cross-sectional p.lane 
of two concentric cylinders is divided into n 
Intervals. 

Outward normal derivative. 

The meridi&n of the free surface in the cross- 
sectional plane is divided into N intervals. 

Gas pressure. 

Liquid static pressure at the height z = 0 . 

The .integral of Q(s) R(s) over the interval. 

A functional of the free surface appearing in the 
linearized Bernoulli equation. 

Radial coordinate. 

Scaled radial coordinate = r/a . 

d/dr . 
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^ 0 ’ 


R(s) 


RO, Rl, R2, ... 


Radii of two concentric right circular cylinders. 

Radius of the equilibrium free surface as a 
fxmction of the arc length along the meridian. 

Normal modes with 0, 1, 2, .., radial nodes in 0, . 


s' Arc length along the meridians of the free surface, 

cylinder wall, and spheroidal bottom in the cross- 
sectional plane: 0 < s ^ S on the free surface. 


subscript s d/ds . 


s, e, Ti 


Surface polar normal coordinates. 




=N 


Midpoints of the N Intervals on the meridian of 
the free surface in the cross-sectional plane i 


t 


Time coordinate. 


t 


3 

Scaled time coordinate = t[(l+jB|)Cf/pa'^]"' 


subscript t 


d/dt . 



Endpoints of the N intervals on the meridian of 
the free surface in the cross-sectional plane. 


T 


Tridiagonal matrix in the discretized time- 
dependent Bernotilli equation. 
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U(z) 


V 


A factor of the velocity potential in the liquid 
contained between two concentric cylinders. 

Fluid velocity. 


V 

W 


Volume of the liquid in the cylinder. 

Matrix in the discrete solution of the Laplace 
equation. 


X(r) 


A factor of the velocity potential in the liquid 
contained between two concentric cylinders . 


y Bessel function of the second kind of order m . 

m 


Vertical coordinate. 


z 


Scaled vertical coordinate = z/a . 


z 


0 


Height of liquid contained between two concentric 
cylinders . 


Z(s) 


Height of the equilibrium surface as a function of 
the radius. 


ZjCr) 


Height of the spheroidal bottom as a function of the 
radius . 


T 

r 


Contact angle. 

Dimensionless growth rate of maximally growing mode. 


AB 


= Aa3^/(da)^/dB) . 
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AB*- 

^critjlO ^crit ’ 

where 

^critjlO 

is the value 


of B . • calculated by the 
crit ’ • 

10-point 

approximation . 

"""on 

~ ^0n,lO ®0n,2O * 




Aoy 

2 2 
= “20 “ ‘*’10 ’ 

4o “10 

are the values 


of - w. calculated by- the 20- and 10-point approxi- 
mations, respectively. 


Ti Displacement normal to the equilibrium surface. 

0 Angie around the cylinder axis. 

00, 01, 02, ... Normal modes with 0, 1, 2, ... angular nodes. 


K 


P 


a 

(j)(r,z,0,t) 

(j>(r,z) 


Capillary constant = pg/P . 

Difference in densities between the liquid and 
gas phases. 

Gas-liquid surface tension. 

Potential function for the fluid velocity. 

A factor in the normal mode expression of velocity 
potential. 


$ 


1 


Vector of values of 0 at N points on the 
meridian on the free surface. 



Vector of values of (j) at M points on the 
meridians on the cylinder wall and bottom. 

Vector = 

Angle in the cross-sectional plane between the 
horizontal and the tangent to the meridian on the 
bottom. 

Angle in the cross-sectional plane between the 
.horizontal and -the tangent to- the meridian- on the 
free surface. 


Frequency of the normal mode. 
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Table A 

2 

1. 03 (B) for various 6 modes; fill height 

radial mode = RO. 

= 0.20; 

-Bond 

ei 

02 

03 

04 

06 

200. 

0.0910 

0.363 

0.817 

1.46 

3.30 

400. 

0.0116 

0.0473 

0.111 

0.208 

0.535 

450. 

0.0027 

0.0121 

0.032 

0.069 

0.226 

500. 

-0.0044 

-0.0161 

-0.031 

-0.043 

-0.020 

550. 

-0.0101 

-0.0392 

-0.083 

-0.134 

-0.223 

600. 

-0.0151 

-0.0585 

-0.126 

-0.210 

-0.392 

700. 

-0.0227 

-0.0889 

-0.194 

-0.330 

-0.657 

800. 

-0.0285 

-0.112 

-0.245 

-0.420 

-0.856 

900. 

-0.0330 

-0.130 - 

-0.285 

-0.491 

-1.01 

1000. 

-0.0366 

-0.144 

-0.317 

-0.547 

-1.14 
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Table 

A 2. o)^(B) 

for variotis 

6 modes; 

fill height 

= 0.30; 



radial mode = RO. 





-Bond 

61 

62 

83 

64 

65 

66 

50. 

0.170 ’ 

0.686 

1.57 

2.86 

4.60 

6.84 

-100. 

0.0328 

.0.139 

0_. 339 

0.666 

1.16 

1.86 

110. 

0.0202 

0.0883 

0.226 

0.464 

0.841 

1.40 

120. 

0.0096 

0.0460 

0.131 

0.295 

0.576 

1.01 

130. 

0.0006 

0.0100 

0.050 

0.151 

0.351 

0.688 

140. 

-0.0072 

-0.0210 

-0.020 

0.027 

0.157 

0.407 

150. 

-0.0140 

-0.0481 

-0.081 

-0.081 

' -0.012 

0.163 

160. 

-0.0200 

-0.0720 

-0.135 

-0.176 

-0.160 

-0.051 

180. 

-0.0302 

-0.113 

-0.225 

-0.336 

-0.410 

-0.411 

200. 

-0.0389 

'-0.147 

-0.301 

-0.468 

-0.612 

-0.701 

210. 

-0.0431 

-0.163 

-0.336 

-0.527 

■ -0.701 

-0.826 

216. 

-0.0462 

-0.175 

-0.358 

-0.562 

-0.750 

-0.892 
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Table A3, w (B) for Various 0 modes; fill height = 
radial mode = Rl. 


-Bond 

00 

02 

04 

06 

50. 

350. 

352. 

359. 

370. 

100. 

153. 

154. . 

158. 

163. 

no. 

135. 

136. 

139. 

144, 

120. 

,119. 

120, 

123. 

127. 

130. 

105, 

106. 

109. 

113. 

140. 

93.5 

94.4 

96,8 

101. 

150. 

82.9 

83.7 

86.1' 

89.8 

160. 

73.3 

74.1 

76.3 

79.8 

180, 

56.1 

56.9 

58.9 

62.1 

200. 

39.6 

40.3 

42.2 

45.1 

210. 

30.1 

30.7 

32.5 

35.2 

216. 

21.1 

21.7 

23.4 

25.8 


0.30 
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Table 

2 

A 4. 03 (B) for various 0 inodes; fill height 

radial mode = RO. 

= 0.40; 

-Bond 

01 

02 

03 

04 

06 

20. 

0.225 

0.936 • 

2.24 

4.28 

11.3 

40. 

0.0384 

0.182 

0.503 

1.10 

3.52 

45. 

0.0161 

0.0920 

0.298 

0.724 

.2.59 

50. 

-0.0024 

0.0172 

0.127 

0.411 

1.81 

55. 

-0.0193 

-0.0498 

-0.025 ■ 

0.133 

1.09 

56. 

-0.0229 

--0.0635 

-0.055 

0.077 

0.938 

57. 

-0.0268 

-0.0784 

-0.088 

0.017 

0.768 

58. 

-0.0329 

-0.100 

-0.134 

-0.069 

0.497 
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Table A 5. 


w^(B) 


for various 8 modes; 


fill height = 0.40; 


radial mode = Rl. 


-Bond 

60 

62 

64 

86 

20. 

132. 

135. 

144. 

158. 

40. 

49.2 

50.8 

55.4 

62.8, 

45. 

38.2 

39.7 

43.8 

50.2 

50. 

28.3 

29.6- 

33.2 

38.9 

55. 

18.1 

19.3 

22.3 

27.0 

56. 

15.8 

16.8 

19.7 

24.2 

57. 

13.0 

14.0 

16.7 

20.7 

• 

CO 

8.07 

8.89 

11.1 

14.3 
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Table A 6 

2 

. w (B) for -variotts 
radial mode = RO. 

6 modes; fill height 

= 0.50: 

-Bond 

91 

02 

03 . 

04 

06 

8. 

0.359 

1.56 

3.93 

7.88 

22.1 

16. 

0.0821 

0.407- 

1.15 

2.51 

7.62 

CO 

H 

0.0455 

0.254 

0.782 

1.78 

5.53 

19. 

0.0281 

0.181 

0.604 

1.42 

4.46 

20. 

0.0082 

0.0993 

0.398 

0.989 

3.09 

20.1 

0.0055 

0.0882. 

0.369 

0.926 

2.88 

20.2 

0.0021 

0.0743 

0.333 . 

0.845 

2.61 

20.25 

-0.0003 . 

0.0645 

0.306 

0.785 

2.40 

20.2759 

-0.0027 

0.0549 

0.280 

0.722 

2.18 
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Table A 7-. O) (B) for various 6 modes;' fill height - 0.50; 
radial mode == El. 


-Bond 

00 

62 

04 

66 

8. 

71.1 

75.8 

89.6 

113. 

16. 

23.3 

'25.8 

32.7 ' • 

44.6 

18. 

16.1 

18.1 

' 24.0 

34.1 

19. 

12.4 

14.2 

19.4 

28.4 

20. 

7.75 

9.24 

13.4 

20.9 

20.1 

7.06 

8.49 

12.5 

19.7 

20.2 

6.17 

7.53 

11.3 

18.0 

20.25 

5.52 

- 6 .. 80 

10.4 

16.8 

20.2759 

4.85 

6.06 

9.41 

15.5 
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Table A 8. w (B) for various 6 modesj fill height = 0.60; 
radial mode = RO. 


-Bond 

01 

02 

03 

04 

06 

3.5 

0.531 

• 2.47 

6.62 

13.8 

39.0 

7.0 

0.152 

0.790 

2.26 

4.81 

13.4 

8.0 

0.0849 

0.487 

1.44 

3.06 

8.16 

8.2 

0.0697 

0.415 

1.24 

2.61 

6.81 

8.4 

0.0483 

0.308 

0.917 

1.88 

4.59 

8.41 

0.0466 

0.299 

0.889 

1.81 

4.40 

8.42 

0.0443 

0.286 

0.847 

1.72 

4.11 
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Table A 

9. w^(B) 

radial 

for various 
mode = Rl, 

0 modes; 

fill height = 0.60; 

-Bond 

00 

02 

04 

06 

3.5 

44.3 

51.5 

73.7 

117. 

7.0 

14.2 

17.8 

28.9 

51.3 

8.0 

7.98 

10.7 

19.0 

35.8 

8.2 

6.39 

8.82 

16.2 

31.0 

8.4 

3.83 

5.70 

11.3 

22.4 

8.41 

3.61 

5.43 

10.8 

21.6 

8.42 

3.30 

5.04 

10.2 

20.6 
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Table A 10. w (B) for various 0 modes; fill height = 0.70; 
radial mode = RO, 


-Bond 

01 

02 

03 

04 

06 

1.5 

0.771 

3.94 

11.1 

23.5 

65.1 

3.0 

0.306 

1.67 

4.79 

9.99 

26.6 

3.2 

0.266 

1.47 

4.20 

8.72 

23.0 

3.4 

0.227 

1.27 

3.63 

7.49 

19.5 

3.6 

0.188 

1.07 

3.06 

6.22 

15.9 

3.8 

0.148 • 

0.855 

2.40 

4.78 

11.8 

3.9 

0.124 

0.719 

1.98 

3.86 

9.26 

3.94 

0.111 

0.647 

1.75 

3.36 

7.93 

3.98 

0.0945 

0.543 

1.42 

2.65 

6,07 



2 

Table A 11, o) (B) for various 0 modes; fill height = 0.70; 
radial mode = Rl. 


-Bond 

00 

02 

04 

86 

1.5 

30,7 

41.5 

78.5 

162. 

3.0 

12.2 

18,1 

38.6 

85.7 

3.2 

10.4 

15.8 

34.5 

77.6 

3.4 

8.67 

13.6 

30.5 

69.2 

3.6 

6.91 

11.2 

26.2 

59.8 

3.8 

4.96 

8,59 

21.1 

48.3 

3,9 

3.75 

6.89 

17,7 

40.7 

3.94 

3.14 

6.02 

15,9 

36.6 

3.98 

2.30 

4.79 

13.3 

30.9 
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